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We investigate Drell-Yan pair production in a QCD inspired model, which takes into account all relevant 
hard processes up to 0(a s ). To address the known shortfalls of such a fixed order calculation we introduce 
phenomenological parton distributions for initial transverse momentum and quark mass, and devise a subtraction 
scheme to avoid double-counting when utilizing the standard longitudinal parton distribution functions. We 
show that we can reproduce Drell-Yan transverse momentum and invariant mass spectra from different proton- 
proton, proton-nucleus and antiproton-nucleus experiments and at different energies without the need for a K 
factor. Fixing our parameters at these spectra, we make predictions for Drell-Yan transverse momentum spectra 
at low hadronic energies, which will be measured for example at PANDA in antiproton-proton collisions. 
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I. INTRODUCTION 

The Drell-Yan (DY) process QJJ] has been studied for the last 
40 years and provides an important tool to access the distribu- 
tion of partons inside the nucleon. While the primary tool for 
exploration of the nucleon structure is deep inelastic scatter- 
ing 10], DY data give complementary insights, since, for ex- 
ample, it directly probes sea-quark distributions y|]. A lot of 
experimental effort is being devoted to measurements of the 
DY process: in antiproton-proton (pp) collisions at PANDA 
(FAIR) [4] and PAX oh, in proton-proton (pp) collisions at 
RHIC fill, J-PARC fflU, IHEP |ll_and JINR JH and in 
pion-nucleon collisions at COMPASS BOU4J] . An overview 
of the experimental situation can be found in IU5I1 . 

Studies of this process lfl6l - [l8ll are generally inspired by 
perturbative QCD (pQCD). The most simple scheme is the 
parton model description, which is a leading order (LO) ap- 
proach (<9(a ( j)). However, it does not fully describe the in- 
teresting observables. While the shape of the invariant mass 
(M) spectra of the DY pair can be reproduced, the absolute 
height can only be accounted for by including an additional 
K factor. Furthermore transverse momentum (pj) spectra are 
not accessible at all [19]. The usual approach to handle the 
latter problem is to fold in a phenomenological Gaussian dis- 
tribution for the parton transverse momentum |20|] , the width 
of which has to be fitted to data. But since these distribu- 
tions are normalized, the absolute size of the cross sections is 
still underestimated [20]. The next logical step is to turn to 
next-to-leading order (NLO, 0(a s )) in the contributing hard 
subprocesses, but this brings about additional problems: the 
calculated pj spectra are divergent for pj — > 0. In fact, they 
are divergent in any fixed order of the strong coupling a s , due 
to large logarithmic corrections In (Mjpj) 11911 . It is possi- 
ble to remove these divergences by an all-order resummation. 
However, since pj is no longer a hard scale at pj — > 0, addi- 
tional non-perturbative (i.e., experimental) input is needed in 
these (and all other pQCD) approaches to describe the region 
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of very small pj 1I21I - I23I1 . Note that the parton model (i.e., 
LO) description is still a very useful starting point, for ex- 
ample for studying spin asymmetries in DY, since there NLO 
corrections appear to be rather small ll24l - l26ll . PANDA, how- 
ever, will allow measurements at hadron cm. energies of a few 
GeV, where non-perturbative effects are expected to become 
more important. This highlights the need to model these ef- 
fects in a phenomenological picture. 

In l27ll we revisited a model, which incorporates phe- 
nomenological transverse momentum distributions for quarks 
and which takes into account the full kinematics in the hard 
LO subprocess, i.e., the usual collinear approximation is over- 
come. It was found that results differed only slightly from 
standard parton model calculations, which underestimate the 
data. In the present paper we improve on our previous work in 
several ways: First, we include quark mass distributions in the 
LO process. We will show that we still underestimate the data. 
This finding triggered a complete calculation of all hard sub- 
processes to 0(a s ) including the full kinematics, which will 
also be presented in the current work. As mentioned above, 
such a calculation would suffer from divergent pj spectra if 
the quarks were massless. However, we will show that the 
phenomenological quark mass distributions we introduced be- 
fore now effectively smear out the divergent behavior. Such 
mass distributions or spectral functions are a well known con- 
cept in nuclear physics, where they are applied to the strongly 
coupled system of nucleons in nuclei, see for example 12811 . 
Thus it is worthwhile to test the same concept in the nucleon, 
which is a strongly coupled system of quarks and gluons. For 
(pr integrated) M spectra the divergent 0(a s ) contributions 
are commonly absorbed into the parton distribution functions. 
Therefore, we introduce a subtraction scheme to prevent dou- 
ble counting of those processes which we consider explicitly. 

This paper is structured as follows: We introduce our nota- 
tion and conventions in Sec. II Al and we give an introduction 
to the general kinematics of DY pair production in Sec. II Bl 
In Sec. [E] we address DY at LO, in particular we describe the 
standard parton model and our extensions of it, namely in- 
cluding intrinsic parton transverse momentum and quark mass 
distributions. We show the results of our extended LO model 
in Sec. ITTEl Sec. [HI] presents our NLO calculation. The ver- 
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tex correction is treated in detail in IIII Al and the calculation 
of gluon bremsstrahlung and gluon Compton scattering is pre- 
sented in Sees. IIII Bl and IIII CI In Sec. IIII PI the influence of 
quark mass distributions is discussed. Furthermore, we dis- 
cuss the treatment of collinear singularities and describe our 
subtraction scheme for the 0(a s ) contributions in detail in 
Sec. IIII El before we shortly comment on the use of initial 
transverse momentum distributions. We show the results of 
our full model in Sec. [TVjand compare with data from differ- 
ent experiments in pp, p-nucleus and p-nucleus reactions. In 
addition we present our predictions for DY pair production at 
PANDA energies. Finally we present our conclusions in Sec. 
[VI The Appendix collects several topics, that are only briefly 
touched in the main text: since we assign different masses 
to the annihilating quarks and antiquarks, we explicitly prove 
gauge invariance in App. [A] Then we investigate the influ- 
ence of different quark masses on the form factors F\ and F% 
in App.|B] and finally we present the details of the phase space 
evaluation in App.|C] 



A. Notation 

In the following we present the conventions and notations 
used throughout this paper: It will turn out to be useful to 
write four-momenta using light-cone coordinates. We employ 
the following convention for general four-vectors a and b: 



a — «o + a z , 
a~ — flo _ o z , 

a± = (a.x, fly) > 
> a 2 = a + a~ — (a±) 2 , 



• b = - (a + b~ + a~b + -2d ± -S x ) . 



(1) 

(2) 

(3) 
(4) 

(5) 



Leptons are treated as massless. We define the target nucleon 
to carry the four-momentum Pi and the beam nucleon to carry 
the four-momentum P2 (see Fig. [TJ. In the hadron center- 
of-mass (cm.) frame we choose the z-axis as the beam line, 
and the beam (target) nucleon moves in the positive (negative) 
direction. Therefore, the nucleon four-momenta read 



Pi = 



Pi 



'Vs 

2 : 

(Vs 



0,0, 



which implies 



PI 



Vs 
2 



trtN— *0 



> Vs 



(6) 
(7) 

(8) 



for the large momentum components of the nucleons. Note 
that in [27] we presented a calculation with vanishing nucleon 
mass Wjyr. However, we want to study our model also at com- 
paratively low cm. energies (e.g. VS ~ 5.5 GeV at PANDA). 



Therefore, in the present work we include the nucleon mass 
since its influence should become significant at these energies. 
We denote the four-momentum of the parton in nucleon 1 (2) 
as Pi (pi)- The on-shell condition in light-cone coordinates 
then reads: 



2 2 + - 
m i = Pi = Pi Pi 



(PiJ 



(9) 



For the virtual photon in Fig. [TJthe maximal q z is derived by re- 
quiring the invariant mass of the undetected remnants to van- 
ish and the photon to move collinearly to the nucleons: 



(Pi + P 2 - qf = X 2 = 

■2Vs > / 9 2 + ( 9z )Lx = o 

S -q 2 



2 VS 



(10) 

(11) 

(12) 



In the literature and in the data presented by many experimen- 
tal groups the Feynman variable 112911 is defined as: 



x F = 



1z 



(<7z)n 



2q z 

VS (<lz)n 

Vs 
2 



(13) 



Note that this approximation for Xp is obviously only valid for 
q 2 <K S . Since we perform studies in the small-S region, our 
definition of the Feynman variable x' F is : 



9z 



(<7z)n 



9z 



2VS 



(14) 



without any approximations. Depending on the experiment 
which we study we will use Xp or x' F according to Eq. (fT3l 
and Eq. <HD . 




FIG. 1. DY production in a nucleon-nucleon collision; Xj and X 2 
denote the nucleon remnants. See main text for details. 
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Pi, mi 




P2-, m 2 r, m r 

FIG. 2. General kinematics of hard subprocesses of DY production. 

B. General kinematics 

In this section we shortly present the kinematic scheme, in 
which our calculations are performed. Note that we show the 
most general form for DY pair production at NLO, of which 
the LO process is a special case. Fig. [2] is a reference for our 
notation. 

For the initial particles we choose four-momenta p\ and p2 
and masses m\ and nt2- For the final state we always choose 
q as the four-momentum of the virtual photon, i.e., the DY 
pair, and M = yfq 2 as its mass. The four-momentum of the 
remaining final state particle we define as r and its mass as 
m r (both at LO). The differential partonic cross section then 
takes the following general form: 

d& = F(p u Pi, q, r) ■ 6 i4 \ P] + p 2 - q - r) dM 2 ±± , 

(15) 

where F contains squared matrix elements, the flux factor and 
constants, and E q (E r ) are the energies of the particles with 
momenta q (r). Note that in the LO case there is only the 
virtual photon in the final state and thus 



F(pup2,q,r) = F{p u pi,q) • 6 (4> {r)2E r dE r 



The Mandelstam variables read 



* = (Pi + Pi) . 
t = (pi- qf , 
u = (pi- qf ■ 



(16) 

(17) 
(18) 
(19) 



In DY measurements the common observables are the invari- 
ant mass M, the absolute transverse momentum pj and the 
longitudinal momentum q z of the lepton pair. Thus we can 
integrate Eq. ( fT~5b over the phase space of r as well as the az- 
imuthal angle of q. Then we obtain 



dfr 



dM 2 df 2^sp c 



F(s, t, M 2 ,m\, m\,m 2 r ) -&{yfs-Eg). 



with the center-of-mass momentum of the incoming states 



y/(s - {mi + m 2 ) 2 )(s - {mi - m 2 ) 2 ) 
Pcm = r-p • (21) 

Now comparing Eqs. ( TT3T > and ( f20b one finds (cf. lfl6ll ) 
E q d& 2 V^Pcm do- 



dM 2 d 3 q ' 
and finally 



n dM 2 dt 



■S({pi+p 2 -q) 2 -m 2 ) (22) 



2 \fsp em {q z ) max do- 



dM 2 dp 2 T dx F 



dM 2 df 



:d((p 1 +p 2 -q) 2 -m 2 ) . (23) 



Thus, for all the relevant subprocesses we actually only have 
to calculate d ^T d? and we will use the relation d23l especially 
for the NLO calculation with two particles (virtual photon + 
quark/gluon) in the final state. 



II. LEADING ORDER DRELL-YAN 

In this section we will present an approach to DY pair pro- 
duction in LO in the hard subprocess, i.e., 0{a { ]). We will start 
with the standard parton model approach and then remedy its 
shortcomings by introducing distributions for transverse mo- 
mentum and mass of the annihilating quarks. Finally we will 
present the results of this LO approach. 



A. Standard parton model 

The LO DY differential cross section in the standard parton 
model reads 0JJ 



dcr L0 = 

I dx x \ dx 2 V q 2 fi(xi , q 2 ) fj(x 2 , q 2 ) d&{x x , x 2 , q 2 ) . 
Jo Jo ; 

(24) 

The sum runs over all quark flavors and antiflavors, q, denotes 
the electric charge of quark flavor i and the functions /; are 
parton distribution functions (PDFs). 

As usual, x\ and X2 are the momentum fractions carried by 
the annihilating partons inside the colliding nucleons, 



Pi = X\P\ 
Pi = x 2 P2 



Note that this implies 



(20) 



Pi 
Pi 



X1 P- 

X 2 PX 



(25) 
(26) 



(27) 
(28) 
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for the large parton momentum components. The small com- 
ponents are 



Pi 



Pi = 



Pi 



(29) 
(30) 



and they are generally neglected in the standard parton model, 
since the partons are assumed to be massless. Consistency 
with the relations ( |25l [26b is then achieved by assuming also 
the nucleon mass to be negligible. 

dcr is the differential cross section of the partonic subpro- 
cess, 



Ana 2 



6(M 2 - q 2 )5 iA \p x +p 2 - q)d 4 qdM 2 , (31) 



and so we find 



dcr 



LO 



Ana 2 



dM 2 dt 9M 2 



S(M 2 -{ Pl +p 2 ) 2 )6(t). 



(32) 



Here q is the four-momentum of the virtual photon, pi, pi are 
the four-momenta of the partons (cf. Fig.[TJ and a « 1/137 is 
the fine-structure constant. 

The factorization into hard (subprocess) and soft (PDFs) 
physics is proven in the collinear case at least for leading twist 
(expansion in 1 /M) in 13011 . 

Note that it becomes immediately clear from Eqs. (l25T l and 
(l26l l that the incoming partons move collinearly with the nu- 
cleons. According to four-momentum conservation in Eq. 
( l3TT l no transverse momentum can therefore be generated for 
the virtual photon (and thus for the DY pair) in the LO pro- 
cess. 

The maximal information about the DY pair that can be 
gained from Eq. (l24T i is double differential, 

dcr L0 v-i 2 ,. ,^^ nal (<7z)max 

= > qffj(xi,M ) f-j(x 2 ,M ) — — 

dM 2 dx F Z-j h ' ja UA ' 9M 2 (P-) 2 E t 



coll 



with 



Xl = 



-(qz)max *F + £coll 



Xl 



coll 



and the energy of the collinear DY pair, 



-■coll 



+ ((<7z)max Xf) 



(33) 

(34) 
(35) 

(36) 



In this section we have presented the standard parton model 
solution for the LO DY cross section. The only quantities in 
this approach not determined by pQCD are the PDFs. These 
have to be obtained by fitting parametrizations to experimental 
data, mainly on deep inelastic scattering (DIS), but also on 
measurements of DY production itself Bill . 



B. Intrinsic transverse momentum 

As already mentioned above, no DY pair transverse mo- 
mentum (pr) is generated in the simple parton model ap- 
proach. Nevertheless, measurements indicate a Gaussian form 
of the pr spectra at not too large pj. This has been studied in 
approaches including initial quark transverse momentum dis- 
tributions, e.g. [20, 32J. In B27I1 we also presented an approach 
incorporating primordial quark transverse momentum to ad- 
dress this issue. However, we found that additional unphysical 
solutions for the momentum fractions x, appear, which have 
to be removed properly. These unphysical solutions are an ar- 
tifact of rewriting the momentum variables using light-cone 
coordinates and they reveal themselves as one of two possible 
solutions of a quadratic equation. The correct solutions can 
always be identified by putting all transverse momenta to zero 
and then by checking whether the well known parton model 
solutions for the x, as given in Eqs. ( 134135b are recovered. In 
l27tl we found that in the transverse momentum dependent ap- 
proach the LO DY differential cross section reads 

do- LO = j~ dxi j- dx 2 J dpi ± J dp 2± 

x ^ Qifiixi . Pl ± . 1 2 )fl(xi, Pi,_ , <7 2 )do-(xi ,pi^,x 2 ,pi ± , q 2 ) . 

(37) 

The functions fj(x, p ± , q 2 ) are now extensions of the standard 
longitudinal PDFs, since they also describe the distribution 
of quark transverse momentum. We will show our ansatz 
for these functions in Sec. HID I Note that we take into ac- 
count the full kinematics in the partonic cross section, i.e. 
do~ lo = d&Lo(xi, pi ± , x 2 , p 2± , q 2 ), however, the quark masses 
are neglected as in the previous section. In this approach the 
transverse momentum (pj = \q±\) of the DY pair is accessi- 
ble, since the annihilating quark and antiquark can have finite 
initial transverse momenta. 

The symbol f represents the requirement that the unphys- 
ical solutions for the x, have to be removed, as mentioned 
above. This is discussed in detail in [27]. Note that in the cur- 
rent paper we take into account a finite nucleon mass. Thus 
the following formulas for the momentum fractions x, are re- 
covered from the corresponding formulas in 12711 by simply 
replacing V5 by P p cf. Eq. ([H}. Then we find for the triple- 
differential hadronic cross section: 



dcr L0 



dM 2 dx 



(P7f 



Fdpi Jo Jo 2 



(<7z)n 



(piJHpiJ 2 



(xit(x 2 ) 2 + (P-) 



F L o((xi)-, pi ± ,(x 2 ) + , p 2± , M 2 ) 
(38) 



with 



^lo((x,)-,/?i ± ,(x 2 ) + ,/ 2± ,M 2 ) = 

X q] fi ((*i)-, K,M 2 ) fl {(xi) + , Pi x , M 2 ) 



Ana 2 
9M 2 



(39) 
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(JCi)- = 



C. Quark masses 



q k ± -q ± 



k±_ ■ q± 



(X2)+ = 



q + k± ■ q± 
— + — + \ 

2 q~ \ 



k±_ -q± 



r 4 



and 



2<7± - k ± > 



1 



E 



2 

max 



q 
q~ 
\q±\ 

k±_ -q± 

(^j_)max 



tJm 2 + pl+ x 2 F (q z ) 

E + X F (q z )max , 
E - x F {q z ) m . ix , 

Pt , 

\Jc±\p T cos((/) ± ) , 

M 2 + p 2 T (l -cos 2 (0 ± )) 



(40) 



(41) 



(42) 
(43) 

(44) 

(45) 
(46) 
(47) 
(48) 

(49) 



In Sec. Ill El we will show and compare the results of this 
approach to the standard parton model and the approach with 
mass distributions as described in Sec. Ill CI 



In Sec. Ill Bl we have extended the standard collinear PDFs 
towards quark distributions which also include quark trans- 
verse momentum. However, we have kept the masses of the 
quarks fixed at zero. Since in light cone coordinates the on- 
shell condition reads m 2 = p 2 — p + p~ - ip±) 2 , this is equiva- 
lent to varying only two of the three, in principal independent 
quark momentum components p + , p~ and p ± . A fully uninte- 
grated parton distribution should depend on all three of these 
components. Therefore, we once more extend the parton dis- 
tributions of Sec. lIIBI bv 



fi(x,p±,q ) 



fi(x,p ± ,m 2 ,q 2 ) . 



(50) 



We will present our ansatz for / in Sec. HID I The differential 
hadronic cross section now becomes 



do-Lo = j~ dxi j- dx 2 j dp 1± j dp 2± J ' dm\ J dm 
x 2 q2 > N Xl ' Pix' ot p q 2 \fl( x 2> PV> m 2> q 2 ) 

i 

X d&(x\,p\ ± , X2, P2 ± ,m\,m\,q 2 ) . 
The partonic cross section is now given by 



(51) 



d<Tio = 



Ana 2 2q 4 — q 2 {m 2 — dm\mi + m 2 ) - (m 2 — m 2 ,) 2 



2 -yjiq 2 - m 2 — m 2 ^ 2 — m 2 m^ 



X 5{M 2 - q 2 ) 6 {4 \pi +p 2 -q) d 4 q dM 2 



(52) 



Note that by assigning different masses m\ and ni2 to the 
quarks, gauge invariance is not preserved at the quark-photon 
vertex. However, the unphysical polarization states of the vir- 
tual photon produced in this manner are projected out at the 
lepton-photon vertex and thus the entire amplitude is indeed 
gauge invariant. See AppendixlAlfor details. 

Now we have to perform the same procedure as described 
in l27ll to remove the unphysical solutions for the longitudi- 
nal momentum fractions Xj. In complete analogy we find (the 
details of this calculation can be found in Appendix lC2b 



dcr 



2_ = d<f> ± -d(k ± ) 2 dm 2 , 

F&Pj JO JO 1 Jo Jo 

(\q± - k*xf + m\ (jq ± +£ ± f + 



dM 2 dx 



dm 2 ^ 



(P-) 2 



(xy) 2 (X2) 2 + {P\) 
X 0(1 -(Jti)_) ((*!)-) ®(1 -(**)+) 0(fe)+) 



Flo((x\)-, pi ± ,m 2 , (x 2 )+, P2 ± ,m\, M 2 ) 



(53) 



with 



Flo((xi)-,pi ± , m\,(x2)+, p2 x ,m\, M 2 ) 



■ Y,q 2 i.h{(M)-,K,mlM 2 ) h{{x2)+,K,mlM 2 ) 

i 

Ana 2 2M 4 - M 2 (m 2 - 6m\m 2 + m 2 ) - (m 2 - m 2 ) 2 



9M 4 



2 -yj(q 2 — in 2 — m 2 ) 2 — m 2 m 



(54) 



2 m 2 
2 



1 



±_ 

2 



k±_ ■ q± 



2 2 
m 2 - m j 



2<? 4 



k± ■ q± m 2 ~ m i 



2\ 



2q+ 



q_ 

q + 



(1 



m 2 + ml ^ 
-Mi-kl- 1 - 2 
4^2 



(55) 
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= — 
r 1 



£+ + m\-m\ 

2 



q 



2q- 



k± ■ q±_ m 2 "M 



and 



2g 



% 1 ^ ^ 



£ = ^M 2 + p 2 + x 2 {q z ) 
q + = E + x F (q z ) m!lx , 

q~ = E - X F (<7z)max , 



2 

max 



|<7j-| = PT , 

k±_ ■ q± - \k ± \p T cos(0j_) , 



iksX 



( m l)max = 



( OT 2) ma x = 



(M 2 +p 2 T )^ 



(56) 

(57) 

(58) 

(59) 
(60) 
(61) 

(62) 
(63) 

(64) 



M 2 + p\{\ - cos 2 ((f> ± )) ' 

2k± ■ q± + q + q~ - J^q + q' {k ± + i<? ± j , (65) 

— 2k ± ■ q ± + m 2 + q + q~ 



Uq + q-m 2 +4q + q-\k ± -^ ± 



(66) 



In Sec. Ill El we show the results of this calculation as well 
as of the approaches in Sees. Ill Al and lHBl 



D. Distributions 

The Bjorken limit and the corresponding infinite momen- 
tum frame, in which the standard parton model is well defined 
and derived from LO pQCD, is an idealization of real experi- 
ments. There the nucleons will always move with some finite 
momentum and thus the partons inside the nucleons will in- 
teract before the collision. These interactions will generate 
momentum components, which are neglected in the (purely 
collinear) standard parton model, namely momentum compo- 
nents perpendicular to the beam line, p 1± ,p2 ± , as well as the 
small light-cone components p\,p 2 - The latter translate to 
non-vanishing quark masses. 

Note that the factorization into hard (subprocess) and soft 
(PDFs) physics is proven in the transverse case at least for par- 
tons with low transverse momentum in li33ll . For the case of 
mass distributions for the quarks we assume this factorization. 



1. Transverse momentum distributions 



the light-cone momentum fraction x,, the transverse momen- 
tum pi ± and the hard scale of the subprocess q 2 . However, the 
general form of these functions is unknown. Known rather 
well are the longitudinal PDFs. Since data of DY pair produc- 
tion are compatible with a Gaussian form of the pj spectrum 
up to a certain pj lf34i[35ll . we assume factorization of the lon- 
gitudinal and the transverse part of /; and make the common 
ansatz |20U3fJ, 13711 



fi(x, p±,q)= fi(x, q ) • fjpj . 



(67) 



Here are the usual longitudinal PDFs and for f x we choose 
a Gaussian form, 



1 



\nD 2 



exp 



(P±Y 

' AD 2 



(68) 



The width parameter D is connected to the average squared 
transverse momentum via 



((pj 2 ) = J dp* ± (pJ 2 fAp±) = 4D 2 



(69) 



and it has to be fitted to the available data. 



Mass distributions 



In Sec. Ill CI we have extended our model by also distribut- 
ing quark masses. This approach is motivated by studies of 
quark correlations and quark spectral functions, see for exam- 
ple IHlllsl]. Note that in such a mass distribution approach 
one effectively parametrizes higher twist effects, i.e. effects 
which are suppressed by inverse powers of the hard scale M. 
These higher twist contributions should become particurlaly 
important in the description of DY pair production in the re- 
gion of small energies (and thus small M), which is one aim 
of our studies. 

For the fully unintegrated parton distributions fi we make 
the ansatz, 



fi(x, p±,m ,q)= fi(x, q ) ■ f±(pj ■ A(p) 



(70) 



Again fi are standard PDFs and f ± are the transverse momen- 
tum distributions of Sec. Ill D 1 1 Since the distribution of lon- 
gitudinal parton momenta is determined by the argument of 
the PDFs x ~ p + , we now allow for a distribution of the re- 
maining degree of freedom, i.e. the small component p~, by 
writing 



A(j>)dp- = 



1 



r(« 2 ) 



N , 



(p~-$) +f 2 (m 2 ) 
with m 2 = p 2 . Rewriting in terms of m 2 yields 



dp- , (71) 



In Sec. Ill Bl we introduced transverse momentum depen- 
dent parton distribution functions ft. They are functions of 
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We choose a non-constant width such that the quark can never 
become heavier than its parent nucleon, 



T(m 2 ) = 



r, 



(73) 



given in several bins of M, xf and p F and for every datapoint 
the average values (M), (xf) and (pr) are given. Since our 
schemes provide Eqs. d38l and d53l l we calculate the quantity 
of Eq. (f75b for every datapoint using these averaged values 
and then perform a simple average in each M 2 bin: 



for < m 2 < m 2 N and f (m 2 ) = otherwise, where T is a free 2E 1 dcr 
parameter. The factor jj normalizes the spectral function such n dx F dp 2 
that 

f dm 2 A(p) = f N dm 2 A(p) = 1 . (74) 
Jo Jo 

where 



do- 



2E 

n VS Jnfi-bm dM 2 dx F dp 2 T 
2E n dcr 



dM 2 



■AM* — — - - «M> , (x F ) , (p T )) 
dM 2 dx F dpi, 



(76) 
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0.1 0.2 0.3 0.4 0.5 

m 2 [GeV 2 ] 



FIG. 3. Spectral function A plotted for different values of the width 
T. Everywhere p + = 0.5 GeV. 



E = ^{M) 2 + {p T ) 2 + {x F ) 2 {{q,) m ^ 2 (77) 

and AM 2 = M 2 w - M 2 in with M max (M min ) the upper (lower) 
limit of the bin. 

We plot the results for the two different approaches of 
Sees. Ill Bl and lH Cl in Fig. [4] The different dashed lines repre- 
sent the massless and the mass distribution approach for dif- 
ferent values of F and they all agree within « 20%. Note, 
however, that with increasing T the calculated cross section 
is slightly enhanced. Everywhere a value of D — 0.5 GeV 
for the transverse momentum dispersion is chosen. With this 
choice of the parameter D the shape of the spectra is described 
rather well. However, in both approaches the absolute size of 
the cross section is underestimated: we have to multiply the 
result of the mass distribution approach for T = 0.5 GeV by 
K = 2 to fit the data (solid line). 



In Fig. [3] we plot A(p) as a function of m 2 for fixed p + = 
0.5 GeV for different values of the width Y. Note that for 
not too small F the region near m 2 — is heavily suppressed 
compared to, for example, T = 0.01 GeV. 



E. Results of the LO calculation 

In this section we present and compare our results for the 
different LO approaches of Sees. Ill Al - IIICI The data are from 
the NuSea Collaboration (E866) lIM and from FNAL- 
E439 1 40]. For the collinear PDFs we used the leading order 
MSTW2008LO68cl parametrizationjil available through 
the LHAPDF library, version 5.8.4 rf42fT 
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1. E866 - pr spectra 

Experiment E866 measured continuum dimuon production 
in pp collisions at S ~ 1500 GeV 2 . The triple-differential 
cross section as given by the E866 collaboration is 



dcr 



2E do- 



d 3 p n y/S dx F dp 2 



(75) 



where an average over the azimuthal angle has been taken and 
where E is the energy of the DY pair, cf. Eq. ( TTTb . The data are 



FIG. 4. p T spectrum obtained at LO from the massless and the mass 
distribution approach with different values of T. Everywhere D = 0.5 
GeV. The solid line is the mass distribution approach for T = 0.5 GeV 
multiplied by a factor K = 2. Data are from E866 binned with 4.2 
GeV < M < 5.2 GeV, -0.05 < x F < 0.15. Only statistical errors are 
shown. 



2. E439 - M spectrum 

Experiment E439 measured dimuon production in pW col- 
lisions at S ~ 750 GeV 2 . The double differential cross sec- 
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tion, 



dcr 



dMdx' F 



(78) 



has been given at a fixed x' F — 0.1. 

As before we begin with Eqs. d38l and d53l > and calculate 
the quantity Eq. ( |78l by integrating over p\ and performing a 
simple transformation from Xf to x' F : 



do- 



dMdx' F jo 



(Pt)1 



dp 2 



do- 



dMdx'dpi 



do- , I M 
— I M, Xf = x' F 1 - — 

dM 2 dx F dp 2 T \ F \ S 



(79) 



Since the experiment was done on tungsten we calculate the 
cross section for pp and pn and average accordingly (74 pro- 
tons and 1 10 neutrons). We compare the results in Fig. [5] Ev- 
erywhere D = 0.5 GeV except for the simple parton model, 
which has no kj distribution. The lowest curve represents the 
indistinguishable results of the standard parton model (Sec. 
Ill Al l and of the (massless) initial kj approach (Sec. Ill Bb . The 
result of the mass distribution approach (Sec. Ill Cb for F = 0.5 
GeV (long dashed) is somewhat larger but still underestimates 
the data: The solid line is the result of this mass distribution 
approach multiplied by a factor K = 1.2 and it fits the data 
very well. 



> 

CD 



0.1 



0.01 



0.001 




r=0.5 GeV, K=1.2 

r=o.5 GeV 
massless 
parton model 



4.5 



5.5 6 6.5 
M [GeV] 



7.5 



FIG. 5. M spectrum obtained from the standard parton model, k T 
approach (massless) and mass distribution approach (D = 0.5 GeV 
for the latter two). The solid line is the result of the mass distribution 
approach multiplied by a factor K = 1.2. Data are from E439 with 
x' F = 0.1. Only statistical errors are shown. 



invariant mass spectra only up to a A" factor and it cannot 
describe transverse momentum (p T ) spectra. The latter issue 
was addressed in the initial kj approach. We found that with 
a suitable choice of an initial kj distribution the DY pj spec- 
tra can be described very well, however, still only up to a K 
factor. The mass distribution approach can improve the pic- 
ture somewhat, but the enhancement of the calculated cross 
sections is too small to describe the data. Still an a priori un- 
determined multiplicative factor K is needed to reproduce the 
measured cross sections. This finding has triggered the NLO 
calculations, which will be presented in Sec. [TTTJ 



III. NEXT-TO-LEADING ORDER DRELL-YAN 

Building on the LO (Oia®)) calculations of Sec. |ll]we here 
present an approach, which incorporates all relevant DY pair 
production processes up to 0(a s ). We will show that by in- 
troducing initial kj as well as quark mass distributions we can 
soften the divergences at low pj of the NLO processes and 
describe pj and M spectra without the need for a K factor. 

In addition to the LO the following processes contribute to 
DY pair production to 0(a s ). First we have the vertex correc- 
tion diagram of Fig. [6] (right). This process alone does con- 
tribute at order a 2 , however, due to identical initial and final 
states it interferes with the LO process of Fig. [6] (left) and the 
interference is of order a s . The same is true for the wave func- 
tion renormalization processes of Fig. [7] Then there is gluon 
bremsstrahlung, where either the quark or the antiquark emits 
a real gluon before annihilating, see Fig. [8] Somewhat dif- 
ferent is gluon Compton scattering since there a gluon and a 
quark/antiquark fuse before or after emitting the virtual pho- 
ton, see Fig. [9] 




FIG. 6. Leading order and vertex correction processes to DY produc- 
tion. Note that only the interference of the two processes contributes 
at NLO. 




FIG. 7. Wave function renormalization processes for DY production 
at NLO. 



F. Conclusion for the LO calculation 

In this section we have presented and compared three dif- 
ferent approaches to DY pair production at LO in the partonic 
subprocess. The standard parton model approach describes 



It is important to note at this point that our initial state 
quarks have different masses m\ and ma when we evaluate the 
matrix elements of the loop correction and the bremsstrahlung 
diagrams, just as in the LO calculation in Sec. Ill CI However, 
we keep the quark mass fixed at the quark-gluon vertex, see 
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Figs. [6] [7] and [8] This guarantees that also in the strong sector 
gauge invariance is preserved, as we show in Appendix[A] For 
the same reason all gluons are treated as massless. 

In the case of gluon Compton scattering we keep the quark 
mass fixed at every vertex, see Fig. [9] for the following rea- 
sons: In principle the final state quark is supposed to be "free" 
and thus one would assign to it a mass ni\ = 0. To preserve 
gauge invariance the quark mass at the gluon vertex must not 
change and thus the exchange quark in the right diagram in 
Fig. [9] would have to be also massless. This, however, im- 
mediately generates an infrared (IR) divergence, as will be 
illustrated in Sec. lIIIDl Therefore, we assign a mass m\ to the 
entire quark line. 




FIG. 8. Gluon bremsstrahlung processes for DY production at NLO. 




FIG. 9. Gluon Compton scattering processes for DY production at 
NLO. 



A. Vertex correction 

1. Formfactors 

The vertex correction process of Fig. [6] together with the 
wave function renormalization of Fig. [7] modifies the bare 
quark-photon vertex: 

y -> = y + sr^(a s ) + o(a 2 j . (80) 

From general principles one now can decompose the correc- 
tion into 

5V = A ■ / + B ■ (p, - p 2 f + C ■ (p l + P2 f, (81) 

where A, B and C are functions of q 2 , nt\ and m 2 . However, 
in DY pair-production the term (pi + p 2 Y does not contribute, 
as we show in detail in Appendix I A II Therefore, from now 
on we neglect this term. 

The Gordon identity l43ll for the case of different masses 
mi and m 2 reads: 

v(p2,m 2 )\ + \u(pi,mi) (82) 

\ mi + nil nil + nil j 

and thus we find: 

io-f v a 

r" = y . (1 + A + (mi + m 2 )B) + — ■ (-(mi + m 2 )B) . 

m\ + m 2 

(83) 

We can now identify the well known form factors F\ and F 2 

Ft = 1 + A + (mi + m 2 )B , (84) 
F 2 = -(mi + m 2 )B . (85) 

The calculation of A and B is tedious, but straightforward. We 
will only need the real parts of A and B, see Eq. (1 1 001 > : 



Re(A) 



a< I t1/ t , \ a I . a + 1 a + <t> + 1 \ 

/ • Re -3 + log(l - v 2 ) - - log(l - a 2 ) + log(l - (a + 4>) 2 )) - - log + log 

\n \ 2 V y 2\a-l a + <p — 1 / 



(b a + d> + 1 
— log 

2 



(86) 



with 



h =- 



1 + - log — 

2 8 1 



1-v 2 



4v 2 + 30 + t + %r I a -l 



2(f> 



2a + < 



log 



a + 1 



+ log 



a + <f> - 1 
a + <A + 1 



+ 2(1 + v 2 + , 



7 3 - 2 log(*) + - log 



(87) 
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h — 

2r 



log(/c) log 



r + 1 - § v 2 -l r+l-S 1 , d> 1 , 

2 ' 1 ~~ 1 - - log 2 (r + 1 - -) + - log 2 (r - 1 - -) 

r-\-t 2 2 2 2 



r- 1 - - 

' 1 2 



+ log — = — log 



2r 



+Li? 



(r+l- 
2r 



Li 2 



(r-l- 
2r 



+ log 2 



1 -v 2 



l - v 2 - : 



+ log 



1 -v 2 



\ 1 - v 2 - 20 



log(tf) 



(88) 



Q-c 1 

Re(B) = -i— Re 
47r mi 



5(1-0;) + ^-/ a-1 ar + 0-l\ t a + 0+1 

— log + log log 

2a + d> \ a+l 6 a + d)+l 4 6 a + d> - 1 



(89) 



and 



2. Soft gluon divergence 



4m 2 



T=(l-V*).|l-=* , 

mi 



1 - 



,2 \ 



'1 ' 



(1-v 2 ) 



' + V z + — , 

4 



(90) 
(91) 

(92) 
(93) 

(94) 
(95) 



Li2 is the Dilogarithm or Spence function. 

We have checked and confirmed that in the limit of equal 
masses m.2 — > mi the well known formula for F\(q 2 ,m 2 ) 
ll44ll45ll is recovered. We note that also for unequal masses the 
ultraviolet (UV) divergences of the loops, displayed in Figs. [6] 
[7] cancel. This is by virtue of the Ward-Takahashi identities, 
which are fulfilled at the quark-gluon vertices, since there by 
construction the mass of the quark does not change. However, 
for m\ + 1112 one finds that gauge invariance is broken at the 
quark-photon vertex (the full amplitude is gauge invariant, see 
Appendix I A Q . The gauge dependence of the quark-photon 
vertex results in a finite renormalization of the charge, which 
cannot be canceled by gauge invariant counterterms. There- 
fore one finds: 



\xmF\(cf,rn\,rn%) + 1 



(96) 



On the other hand, q 2 - M 2 sets the hard scale in DY pair pro- 
duction and thus our model should only be applied at reason- 
ably large q 2 . A sensible lower limit would be q 2 > 1 GeV 2 . 
Thus we probe F t far away from q 2 = and we show in Ap- 
pendix|B]that for these physically interesting q 2 the influence 
of the renormalized charge is negligible. 



To obtain Eqs. ( |86ll95l l we have assigned to the gluon a mass 
A which serves as an IR regulator in the loop integral. Ac- 
cording to the theorems by Kinoshita-Poggio-Quinn J4^-SH 
and Kinoshita-Lee-Nauenberg [46, 50:] the (IR) divergence of 
the loop integral cancels against a similar divergence of the 
bremsstrahlung processes in Fig. [8] Therefore we will also 
introduce the same gluon mass A in the calculation of the 
bremsstrahlung in Sec. IIII Bl and we will show in Sec. IIV Al 
that the two divergences actually cancel numerically, as they 
should. 



3. Interference cross section 

The LO partonic cross section can be written as (M 2 = q 2 ) 
da"LO 1 nor 2 



dM 2 dt 3 6M 3 pc 



■T h0 -6(s-M z )6(t-mi), (97) 



where I is the color factor and 



'LO 



g,v + ^ I • Tr [q» 2 - m 2 yf^> x + mi )/] . (98) 



From this one easily finds that the cross section of the inter- 
ference of the LO process and the vertex correction can be 
written as 



do- 



ve 



1 no 2 



c\M 2 c\t 3 6M 3 Pc 

with 

Tvc = -gftv + — 2~ 



T vc -S(s-M z )6(t-m 2 ) . 



(99) 



x (2Re(A) • Tr [(f 2 - m 2 )y»(f> 1 + m x )y v ] 

2 Re(fi) ■ Tr [(f 2 - m 2 )/(^ + mi)] • (p 1 - p 2 ) v ) . 



(100) 



Note that dtryc depends on the gluon mass A of Eq. 
and so does the cross section for gluon bremsstrahlung of 
Sec. IIII Bl Only the sum of the two cross sections is a physi- 
cally meaningful quantity and thus we will only plot the sum 
of the two in our results in Sec.lIVI 
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Kinematics 



D. Influence of quark mass distributions on DY p T spectra 



Since the vertex correction shares initial and final states 
with the LO process, the hadronic cross section of the vertex 
correction is calculated exactly as described for the LO mass 
distribution case in Sec. Ill CI 



B. Gluon bremsstrahlung 

In the case of gluon bremsstrahlung we assign the same fic- 
titious mass A to the gluon as for the vertex correction process. 
This ensures the cancellation of the soft gluon divergences, as 
we will show in Sec. [TV] Then the partonic cross section be- 
comes (E g > A) 



4 a 2 a s 



do-B 



dM 2 df 9 48M 2 sp 2 cn 



®(E g ) . 



(101) 



Here ^ is the color factor and Tg is given by 



with 



Ti[(f 2 -m 2 )S a ^ 1 +m 1 )S v a ] (102) 



fx -Q + mi 



i-ti 



+ mi 



, < , r tr- (103) 

(pi - q) - '«2 (pi ~ 9) ~ m l 

The calculation of the hadronic cross section basically fol- 
lows along the same lines as for the LO case in Sec. Ill CI Once 
again one has to remove unphysical solutions for the momen- 
tum fractions x,, however, the calculation of the phase space 
is more subtle. The details of this calculation are given in ap- 
pendix |UT] 



C. Gluon Compton scattering 



As already mentioned in the introduction, massless pQCD 
calculations of DY pair production at NLO produce diver- 
gent pr spectra [19]. The origin of these IR divergences are 
twofold: first there is soft gluon emission (bremsstrahlung). 
This type of soft divergence, however, is not problematic, 
since it exactly cancels against a divergence in the virtual 
processes (vertex correction), cf. Sec. IIII A 21 Second there 
is emission (bremsstrahlung) or capture (Compton scattering) 
of a gluon by a massless participant quark (also called mass 
or collinear singularity): the M-channel exchange quarks in 
Figs. [8] and [9] can become onshell at pj — and thus for 
m\ — m.2 = the propagators in Eqs. (11031106b produce a 
non-integrable (in p T ) singularity at pj = 0. To address this 
problem was one reason for introducing mass distributions for 
the participating quarks. This procedure aims at smearing out 
the divergence and so making the pj spectra integrable. 

We will illustrate this procedure on the example of the 
gluon Compton scattering process. In Fig. [10] we compare 
pj spectra produced by gluon Compton scattering in two dif- 
ferent schemes: one calculation with massless quarks and a 
calculation which includes quark mass distributions. In both 
cases the quark's initial transverse momentum is set to zero. It 
is seen that now indeed a transverse momentum of the dilep- 
tons is generated. However, its magnitude is significantly un- 
derestimated. One can also see clearly that the rise for pj — > 
of the calculation with mass distributions is slower than for the 
calculation with massless quarks. This is a consequence of the 
effective cut-off, which is introduced by distributing the quark 
masses. We find that the divergence in pj is softened enough 
to make the pj spectra integrable. 



For gluon Compton scattering we choose for the initial 
quark/antiquark to have four- momentum p \ and for the gluon 
to have four-momentum p2, thus ;«2 = since the gluon is 
real. For the outgoing quark/antiquark we then have m r = 
The partonic cross section then reads (E r > m r ): 

d&c 1 a 2 a s 

T c ■ ®(E r ) , 



dM 2 df 6 12M 2 (s - m 2 ) 2 
where 4 is the color factor and Tc is given by 



?c = - 
with 



Tr 



+ pr,~(t + m 1 )S 



/'<>/- 



• OTl 



(104) 

(105) 
■f . (106) 



(pi + pi) 1 - m\ (p\ - q) 2 - m\ 

The calculation of the hadronic cross section is similar to 
the case of gluon bremsstrahlung, however, the inherent asym- 
metry of the initial state (massive quark hits massless gluon) 
requires additional care. The details can be found in appendix 
IC31 
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FIG. 10. p T spectrum of gluon Compton scattering obtained from 
massless and mass distribution approach with initially collinear 
quarks. The PDFs are the MSTW2008LO68cl set. Data are from 
E866 binned with 4.2 GeV < M < 5.2 GeV, -0.05 < x F < 0.15. 
Only statistical errors are shown. 
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E. Collinear (mass) singularities and parton distribution 
functions 

As we have just shown, we have regularized the collinear 
singularities of the NLO processes by introducing quark mass 
distributions. However, for the calculation of the cross sec- 
tions we would like to use PDFs as supplied in the litera- 
ture. But exactly those collinear singularities, that we have 
just regularized, are commonly absorbed into the definition of 
the standard PDFs. Thus, to avoid double-counting, in this 
section we present a subtraction scheme, that leaves us with a 
consistent cross section to the order of the hard processes we 
are considering. 

To set the stage we first review briefly the introduction of 
the renormalized PDFs into pQCD. For the DY process this 
concerns the calculation of M spectra, because the spectra dif- 
ferential in pr are not accessible by pQCD. Since we are in- 
terested in the description of fully differential DY spectra by 
our model, we have to modify the way towards the standard 
renormalized PDFs. This is outlined in a second subsection. 



1. Collinear singularities in pQCD 

If Bjorken-scaling were not violated the PDFs found in 
deep inelastic scattering were functions of the momentum 
fraction x only. However, it is well known that the interac- 
tions among the quarks and gluons induce scaling violations 
via processes like gluon bremsstrahlung and gluon quark- 
antiquark production. To calculate the contributions of these 
processes to the longitudinal PDFs one has to integrate over 
the transverse momentum of the emitted particle (quark or 
gluon) and finds, that they suffer from collinear (or mass) sin- 
gularities, i.e. they are singular because the quarks are treated 
as massless. The divergences appear at the boundaries of the 
transverse momentum integrals (which is why they are called 
collinear divergences). Thus, one can regulate these diver- 
gences by introducing a regulating cut-off rf in the transverse 
momentum integral. In this scheme one can define renormal- 
ized longitudinal quark PDFs by absorbing the collinear sin- 
gularities and the (non-measurable, scaling) bare quark and 
gluon PDFs, jf(x) and g°(x), into one function lUlll : 



fi(x,f) 



TV) 



+ g°(y) 



(107) 



with the hard scale fj 2 . The coefficient functions of the diver- 
gent logarithms are the splitting functions P qq and P qg , which 
are given below. The functions C q and C~ contain possible 
finite contributions of the scaling violating processes and the 
superscript S reminds us of the fact, that these finite contri- 
butions depend on the chosen renormalization scheme, since 
only the divergent contributions actually have to be absorbed 



into the renormalized PDFs. The splitting functions found in 
deep inelastic scattering to order a s are given by l5lll 

v-2 

X) 



1 +. 



1 r 



(l-x) 

(l-x) 2 ] , 



(108) 



(109) 



where | and i are color factors. The plus prescription reads: 



Jo (1-*)+ Jo 



d.v 



/( *)-/(!) 
l-x 



(110) 



where / is a smooth function on [0, 1 ] . 

Remarkably one finds the same splitting functions when 
calculating order a s corrections to DY pair production 15 ill : 
P qq collects all the contributions from the processes with 
qq in the initial state, i.e. the vertex correction and gluon 
bremsstrahlung processes of Figs. [6] [7] and [8] P qg contains the 
contributions from the gluon Compton scattering processes in 
Fig-E] Then the partonic cross section for DY pair production 
to order a s integrated over the transverse momentum of the 
DY pair can be written schematically as ll5~lll 



dcr Ana 



AM 2 9M 4 
for a quark of flavor i and with 

M 2 M 2 



z = 



XIX2S 



X\X2 



(HI) 



(112) 



where y[S ( yS) is the partonic (hadronic) cm. energy and 
X\,X2 the momentum fractions of the quarks. The <5-function 
in d 1 1 II ) gives just the leading-order contribution and the func- 
tions T qq and T qg give the contributions with initial states 
consisting of quark-antiquark (vertex correction and gluon 
bremsstrahlung) and quark-gluon (gluon Compton scattering), 
respectively: 



, M 2y 

T qq (z) = 2P qq {z)log\-^\-\ O.), 



T qg (z) = P qg {z)log 



M 2 



(113) 



(114) 



where again the cut-off rj 2 was introduced to regulate the trans- 
verse momentum integration and where the functions C q and 
C g are again renormalization scheme dependent, finite con- 
tributions. In principle one could obtain the hadronic cross 
section by folding the partonic cross section ( 111 II ) with the 
bare parton distributions and summing over all quark flavors: 

dcr 
dM 2 

Ana 2 , f 1 , , t 

dxjdx2 &(x\X2 — t) 



9M 4 



X\X2 



{(jf(xi)jf(x2) + <jf <-> /f))[<5(i 
(gVo (f?(x 2 ) + jf(x 2 )) + (/ » jf,jfi) ^r qg (z) 



■z) + -r q , 



(-) 



(115) 
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This cross section cannot be evaluated straigtforwardly, since 
neither the bare parton distributions are available, nor are T q q 
and T qg well defined, since they depend on the arbitrary cut- 
off r] 1 . However, we note the following relation for a general 
function P: 



f d Xl dx 2 — (JC1JC2 - T)ff{x l )fi. , {x 2 )P 
Jo X\X2 \.V|.V 2 

Jt Xl Jr. X 2 \X\X2j 

*rW~K)j>(?)* 

t J dxidx 2 dz5(l - z)6(x 1 x 2 z - T)ff(xi) 



?<y) 



(116) 



In addition one finds for the product of two renormalized 
PDFs of the type of Eq. JTOTh 



=ff(x l )f?(x 2 ) 



+ 8°W 



M?M?H(?) 



M 2 



log|^ +C q 



s I x ± 

y 



+ g°(y) 



+o(«J) • 



(117) 



Comparing the last two equations, one finds, that one can ex- 
press the hadronic cross section d 1 15b in terms of the renor- 
malized PDFs ( TUFTY One obtains [51] 



dcr 
dM 2 



4na 2 sr- 
"qIF T Zj 



dx\dx 2 dz6 {x\x 2 z - t) 



9M 4 

x {(/;■(*!, M 2 )/K* 2 ,M 2 ) + (/; <-» / 7 )) 
+ [*(*!, M 2 ) (/ife, M 2 ) + / 7 (* 2 , M 2 )) 



^(l-z) + ^C'(z) 



+ 0?" ft, ft] ^C s g (z) 



(118) 



which is correct to 0(a s ). The collinear divergences and the 
bare PDFs have been absorbed into the renormalized PDFs. 
Again there remain scheme-dependent finite contributions C„, 
C e . 



2. Collinear singularities in our model 

The introduction of the renormalized PDFs gives rise to the 
famous DGLAP evolution equation which describe success- 
fully the scaling violations 15111 . Clearly for our model we 
want to inherit this fundamental QCD property. Also from a 
pragmatic point of view we want to use the standard (renor- 
malized) PDFs from the literature. On the other hand, in 
the pQCD approach to the DY process the reshuffling of the 
collinear singularities into the renormalized PDFs is only pos- 
sible for the pt integrated M spectrum. In our model we also 
want to describe the DY spectra differential in pj. Therefore, 
we cannot follow the steps outlined in the previous subsection. 
However, we have an explicit regularization of the collinear 
singularities. This allows to make contact between the bare 
and the renormalized PDFs by a kind of backward engineer- 
ing, which we will describe next. 

Since we explicitly take into account all 0(a s ) processes, 
cf. Sees. IIII A11III Cl we now have to ensure that all the 
0(a s ) contributions to the cross section, that were absorbed 
into the renormalized PDFs are subtracted. Otherwise we 
would double-count the 0(a s ) contributions. Schematically 
the hadronic cross section for DY pair production to next-to- 
leading order in a s can be written as 



dcr 



derm f?-f?+ddyc+Bjf-j? 

LO(a?)=0(l) 0{a s ) 

+ d& c g°-(jf + jf) , (119) 

0(a s ) 

with the bare PDFs jf, g°. We note that /; ■ fi = jfjf + 0(a s ) 
and g-(f i+ fi) = g° ■ iff + Jf) + 0(a s ), and 'so' to 0(a s ) 
nothing changes if we replace the bare PDFs multiplying the 
NLO partonic cross sections: 



dcr 



d^o ^./;°-d^(_B ./;•./; 

0(a?)=0(l) 0(a,) 



+ d(Tc g-(fi+fi) 

0(a s ) 



+ 0(a l s ) . (120) 



However, if we were to do the same with the LO term, we 
would get additional 0(a s ) contributions, cf. Eq. dl 171) . How 
do we subtract these contributions? After all, we cannot cal- 
culate the integrals in Eq. dl 17b . since we do not know tj or 
the bare PDFs. 

In our model we set out to calculate exactly those transverse 
momentum spectra that were integrated in the derivation of 
the renormalized PDFs in the pQCD case above. To accom- 
plish this, we have introduced quark mass distributions to han- 
dle the collinear divergences that enter the quark PDF in Eq. 
(1107| i. The important difference between the pQCD approach 
and our model is, therefore, the following: in pQCD the regu- 
lating cut-off rj 1 is completely arbitrary and physical results 
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can only be obtained by "absorbing" this arbitrariness into 
the renormalized PDFs. All the finite and scheme-depedent 
contributions can be calculated analytically. In our model we 
know the regulator 77 in ( 1107b . it is nothing but our quark mass 
m (or better m 2 ). However, we do not know the finite contri- 
butions Cq „ in our model. To estimate them, we introduce two 
new parameters K q and K g , so that the functions T q q and T qg 
become: 



T?M=2P q J Z )\og 



( M 2 



v K 2 m 2 



r q m g (z) = p qg (z)io g 



( M 2 



(121) 
(122) 



Then we can rewrite the renormalized PDFs in Eq. ( II 171 ) as 



fi(x u M 2 )fi(x 2 ,M 2 ) 



=jf(x l )fi°(X 2 ) 



+fi°(x l ) 



In 
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Solving for the product of the bare PDFs gives 



(123) 



jf(xi)jf(x 2 ) 
=fi(xuM 2 )f 1 (x 2 ,M 2 ) 
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P qg \ 
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^K 2 m 2 ) 



+0(a 2 s ) . 



(124) 



On the right hand side we can replace all the bare PDFs by 
their renormalized version, since all additional corrections in- 



troduced by this procedure are 0(a 2 ): 
ff{^{x 2 ) 



=fi(x 1 ,M 2 )f i (x 2 ,M 2 ) 



-ftxuM 2 )^ 
2n 



M 2 ) 



«| — l lo 8 



< M 2 ^ 



2 2 
\ K q m J 



+ g(y,M 2 ) 



M 2 



-/;(*2,M 2 ) 

+o(«5) . 

We define 



yK 2 s m 2 j 



fi(y,M 2 ) 



p. l9 \j\lw 



M 1 



yK 2 m 2 j 



+ g(y,M z ) 



< M 2 ^ 



(125) 



fif°(x, M 2 ) = g £ dy i^Cy, M 2 )P W ) log 

f dy ig(y, 
2tt J x y 



g suh (x, M 2 ) 



and so we find 



M 2 )P qg \-_\\o % 



' M 2 ^ 
K K 2 m 2 j 
( M 2 ^ 



K K 2 g m 2 j 



(126) 



ff(M)jf(x 2 ) 
=f i (x 1 ,M 2 )fi(x 2 ,M 2 ) 

-fi( Xl , M 2 )f{ ub (x 2 , M 2 ) - fi( Xi , MV ub (x 2 , M 2 ) 
-/](x 2 , M 2 )/; sub (^!, M 2 ) - Mx 2 , M 2 )g™\ Xl , M 2 ) 
+O(o?) • (127) 

In this scheme the hadronic cross section dl20t becomes 

- do-, () (fi . ./;;- + ■ ./;) + d<x vc+ B ft ■ fi 
-dfr L0 (fi ■ g™ h + fi ■ g suh ) + d& c g- (fi + fi)] 

+ 0(a 2 ) . (128) 

From now on, we label the different contributions to the cross 
section as sketched in the following: 



dcr L0 = 



(129) 



dVqq = j ^ ^ d ^VC + B fi ■ fi ~ drr LO (fi ■ ff^ + j?*> ■ fi) , 

(130) 

J ^e*d& c g ■ (fi + fid 

-dfr LO (fi-g™ h + frg™ h ). (131) 



dcr 
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In Eqs. ( 1130b and (1131b we subtract now precisely those 0(a s ) 
contributions which were absorbed before into the renormal- 
ized PDFs. Thus, by this procedure we have produced in our 
model a consistent cross section to 0(a s ). Indeed, we have 
checked explicitly that the collinear divergences which appear 
in d<xvc+B, if the quark masses are sent to zero, cancel the cor- 
responding divergence of ff uh as given in Eq. (1126b . The same 
is true for dcr c and g sub . Thus the expressions (1130b and ( 1131b 
remain finite for vanishing quark masses. 

For the calculation of our results in Sec. lIVI we will use Eqs. 
(1129b - (1131b . The quantities which enter are, on the one hand, 
the standard PDFs which we can take from the literature, and, 
on the other hand, the parameters K q and K g which appear in 
Eq. (1126b . We recall that these parameters introduced in Eqs. 
(1121b and (1122b correspond to finite, i.e. infrared safe, con- 
tributions which appear, e.g., as C qt , in Eqs. (11 13b and (II 14b . 
We will vary the parameters K q and K g around natural values 
(i.e. around 1), to estimate these finite contributions. 

F. Initial transverse momentum distributions 

Even in the mass distribution approach at NLO we find that 
the pr data are heavily underestimated. Therefore, we have 
also introduced parton initial transverse momentum distribu- 
tions for the NLO processes, just as we did at LO in Sec. Ill B I 
Taking into account these distributions we obtain a good de- 
scription of measured cross sections without K factors, as we 
will show in Sec.lIVI 



groups, for example ll4ll l52Tl . and they are not unique. We 
will show results obtained with different leading order PDFs 
to get an estimate for the theoretical uncertainty. Once again 
we employ PDFs available through the LHAPDF library, ver- 
sion 5.8.4 JH. 

At this point we note two things: first, the divergence of the 
NLO processes near pj — can also be cured by choosing a 
finite and fixed quark mass. However it is much more natural 
to assume a broad mass distribution, since the nucleon is com- 
posed of strongly interacting partons; we will show below, 
that we can describe the experimental data quantitatively well 
assuming broad quark mass distributions. Second, in princi- 
ple we could have chosen different transverse momentum and 
mass distributions for valence quarks, sea quarks and gluons 
(only transverse momentum in this case). We have not done 
so for two reasons: on the one hand this would have forced 
us to introduce additional parameters, while one is of course 
always inclined to keep the number of parameters as small as 
possible. On the other hand, as one will see in the results, the 
data do not require additional modelling. Thus we always use 
the same distributions from Sec. lIIDI for all quarks and gluons 
(again no mass distributions for the latter). 

A. E866 - p T spectrum 

In this section we present the results of our full model. 
The data are from E866 IT341 l35tl for pp collisions at S = 
1500 GeV 2 . 



IV. RESULTS 

In this section we present our combined results of the LO 
and NLO calculations of Sees. [IT] and [Til] Since one of our 
main goals is to gain access to the DY pj spectra, we compare 
our model with data on pj spectra from different experiments. 
Most of these experiments were done in proton-nucleus col- 
lisions. However, the recent E866 experiment 1134 13511 mea- 
sured pp collisions and we will begin to fix the parameters 
of our model at the data for this more elementary reaction in 
Sec. IIV Al Using the fixed parameters we compare our results 
to data for different pp and pA experiments in Sees. IIV BlIIV Fl 
In Sec. lIVGl we study our results for an antiproton-nucleus ex- 
periment (pN) (E537), which is useful, since one of our aims 
is to make predictions for DY pair production at PANDA, 
which will measure antiproton-proton (pp) collisions. These 
predictions can be found in Sec. IIV HI 

Our model has effectively four parameters: The width of the 
initial parton transverse momentum distribution, represented 
by D, the width F of the quark mass distribution (spectral 
function) and the two subtraction constants K q and K g which 
estimate the finite contributions which were absorbed into the 
PDFs. We will vary K q and K g in a natural range | . . . 2. The 
gluon mass A is not really a parameter but a numerical neces- 
sity and we will explore its influence on our calculations in 
Sec. IIV Al In addition our model utilizes standard parton dis- 
tribution functions. These are calculated by several different 
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FIG. 11. pr spectrum obtained from our full model with different 
values of D. Everywhere T = 0.2 GeV, A = 5 MeV, K q = 1 and 
K g = 2. The PDFs are the MSTW2008LO68cl set. Data are from 
E866 binned with 4.2 GeV < M < 5.2 GeV, -0.05 < x F < 0.15. 
Only statistical errors are shown. 

First we want to stress, that the calculation in our full model 
reproduces measured DY pj spectra without the need for a K 
factor, see for example Fig. [14] In order to better understand 
the parameter dependence of our model, in the following we 
explore the parameter space. The details of how we obtain the 
presented cross sections were given in Sec. Ill E 1 1 First we 
show several plots, in which we vary only one parameter and 
keep the others fixed. 
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In Fig.QT|we plot the results of our full NLO model for dif- 
ferent D. As one can see for D around 0.45 GeV, which corre- 
sponds to an average squared initial quark transverse momen- 
tum (£ 2 ) = (0.9) 2 GeV 2 , the data are reproduced quite well. 
Obviously, D determines the shape of the pj spectra. 
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FIG. 12. p T spectrum obtained from our full model with different 
values of T. Everywhere D = 0.45 GeV, A = 5 MeV, K q = 1 and 
Kg = 2. The PDFs are the MSTW2008LO68cl set. Note that the 
curves for T > 0.2 GeV are on top of each other. Data are from E866 
binned with 4.2 GeV < M < 5.2 GeV, -0.05 < x F < 0.15. Only 
statistical errors are shown. 

In Fig. [12] we show results for different V. The results for 
several values of F all agree very well with each other and at 
the same time reproduce the data quite well. Thus at E866 
energies our model appears to be rather insensitive to changes 
of T over a wide range. 
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FIG. 13. pr spectrum obtained from our full model with different 
values of A. Everywhere D = 0.45 GeV, T = 0.2 GeV, K q = 1 and 
K g = 2. The PDFs are the MSTW2008LO68cl set. Note that the 
curves for A < 5 MeV are on top of each other. Data are from E866 
binned with 4.2 GeV < M < 5.2 GeV, -0.05 < x F < 0.15. Only 
statistical errors are shown. 

Remember that the gluon mass A is introduced to regulate 
divergences that occur when the gluons become very soft. In 
the limit of A — > these divergences of the bremsstrahlung 



and the vertex correction processes should exactly cancel. 
Therefore, it is sensible to choose A as small as numerically 
feasible. Results for different choices of the gluon mass A are 
shown in Fig. [13] While the results for A = 100, 250 MeV are 
still visibly larger than the results for 5 and 0.5 MeV, the latter 
two agree very well with each other. Thus the calculated cross 
section appears to converge in the A — > limit. Therefore, we 
chose A = 5 MeV for the calculation of all the following re- 
sults. 
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FIG. 14. p T spectrum obtained from our full model with different 
PDF sets. Everywhere D = 0.45 GeV, T = 0.2 GeV, A = 5 MeV, K q = 
1 and Kg = 2. Data are from E866 binned with 4.2 GeV < M < 5.2 
GeV, -0.05 < x F < 0.15. Only statistical errors are shown. 



We show an example of the influence of different par- 
ton distribution functions in Fig. [14] The results with the 
MSTW2008LO68cl JH and GJRO8I0 [52] sets agree quite 
well with each other, with only small deviations. This is an 
illustration of the uncertainties induced by the different inte- 
grated parton distributions. 
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FIG. 15. pt spectrum obtained from our full model with differ- 
ent values of the subtraction parameter K q . Everywhere D = 0.45 



GeV, T = 0.2 GeV, A = 5 MeV and k. 



2. The PDFs are the 



MSTW2008LO68cl set. Data are from E866 binned with 4.2 GeV 
< M < 5.2 GeV, -0.05 < x F < 0.15. Only statistical errors are 
shown. 
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FIG. 16. pr spectrum obtained from our full model with differ- 
ent values of the subtraction parameter n g . Everywhere D = 0.45 
GeV, T = 0.2 GeV, A = 5 MeV and K q = 1. The PDFs are the 
MSTW2008LO68cl set. Data are from E866 binned with 4.2 GeV 
< M < 5.2 GeV, -0.05 < x F < 0.15. Only statistical errors are 
shown. 



To determine the subtraction parameters K q and K g we ex- 
plore the dependence of the cross section on these two pa- 
rameters in Figs. n"5land[T6l In the range of natural choices 
(K q , Kg - j . . . 2) we find, that with K q - 1 and K g — 2 the data 
are desribed rather well. Although Fig. [15] would indicate a 
better fit for smaller K q , we stick to K q = 1, since for this value 
the slope of the M spectra fits also very well, as we show for 
example in Sec. II V Dl 



the cross section. The contribution of the qg correction is rel- 
atively small and negative for not too large pj. 

Fig.[T"8"1 shows a comparison of the results of our model with 
the results of a PYTHIA event generator calculation. For this 
specific plot PYTHIA version 6.225 JH with CTEQ5L PDFs 
15411 is used and the results are scaled up by a factor K — 2. 
The comparison with the experimental data suggests a good fit 
of the shape of the spectrum for a value of (k^) = (0.8 GeV) 2 
in PYTHIA (internal parameter PARP(91) = 0.8). For out- 
value of D we get {k 2 T ) = (2D) 2 = (0.9 GeV) 2 . Although the 
PYTHIA parameter is obviously intended to have the same 
meaning as our definition for the average initial transverse 
momentum, the complex internal treatment of the interaction 
in the PYTHIA code leads to some numerical mismatch with 
our implementation. 
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FIG. 17. p T spectrum obtained from our model decomposed into the 
different contributions as described at the end of Sec. lIIIEl Note that 
we plot the negative quark-gluon contribution. Everywhere D = 0.45 
GeV, T = 0.2 GeV, A = 5 MeV, k„ = 1 and Kg = 2. The PDFs are the 
MSTW2008LO68cl set. Data are from E866 binned with 4.2 GeV 
< M < 5.2 GeV, -0.05 < x F < 0.15. Only statistical errors are 
shown. 



In Fig. [T7] the cross sections of the different contributions 
to the full result are plotted individually. The definition of 
the LO, quark-antiquark (qq) and quark-gluon (qg) contribu- 
tions is given at the end of Sec. IIIIEI Note that the sum of 
the LO contribution and the qq corrections make up most of 



FIG. 18. pr spectrum comparison of PYTHIA results with our full 
model. For the PYTHIA calculations version 6.225 with CTEQ5L 
PDFs were used. The PYTHIA result is plotted for two different 
values of the average initial k T and multiplied by a factor K = 2. Our 
model was calculated with D = 0.45 GeV, T = 0.2 GeV, A = 5 MeV, 
K q = 1, K g = 2 and MSTW2008LO68cl PDF set. Data are from E866 
binned with 4.2 GeV < M < 5.2 GeV, -0.05 < x F < 0.15. Only 
statistical errors are shown. 



For comparison we show our results for a different M bin in 
Fig. [19] With the parameters determined above our full model 
reproduces the measured spectrum very well. In contrast to 
the LO approach of Sec. [II] we do not need a K factor to de- 
scribe the absolute height of the spectrum. At the same time 
the width D of the intrinsic (non-perturbative) kj distributions 
changes only little (D = 0.5 GeV vs. D - 0.45 GeV) when 
passing from the LO to the NLO calculation. Note, however, 
that Fig.[l7]indicates, that at least the results for the qq correc- 
tions deviate from a Gaussian behavior at pj > 3 GeV. Thus 
the contributions of some of the hard NLO processes are only 
significant in the high (perturbative) pj region. Therefore, to 
describe the low pj regime at NLO one still requires almost 
the same non-perturbative input for the intrinsic parton kj, i.e., 
only little transverse momentum is generated dynamically and 
the width parameter D does not change considerably. 



18 



1 

> 
to 
ED 

Q. 



"D 

o 

"D 
LU 



0.1 



0.01 




0.5 1 1.5 2 2.5 3 3.5 4 
Pt [GeV] 



FIG. 19. pr spectrum obtained from our full model. Everywhere 
D = 0.45 GeV, T = 0.2 GeV, A = 5 MeV, K q = 1 and k s = 2. The 
PDFs are the MSTW2008LO68cl set. Data are from E866 binned 
with 7.2 GeV < M < 8.7 GeV, -0.05 < x F < 0.15. Only statistical 
errors are shown. 



E 2 = M 2 + (p T ) 2 m . dx + 4(?z)max 
2n2 



(PtX 



(S +M Z - M Z R ) 
45 

(S + M 2 - M 2 ) 2 
4S 



- M 2 - 4fe)^ 



(137) 
(138) 

(139) 



B. E866 - M spectrum 



where M 2 R is the minimal invariant mass of the undetected 
remnants. For pp and pn collisions we choose a value of 
Mr - 2niN and for pp a value of Mr - 1.1 GeV, which can 
be interpreted as two times a diquark mass. Note that at cm. 
energies of Vs" * 15.3 GeV (E537), Vs" ~ 27.4 GeV (E288, 
E439) and Vs -38.8 GeV (E605, E772, E866) we are not 
really sensitive to these values if they stay at or below a few 
GeV. 

In Fig. |20]we compare our result for the double differential 
cross section to the data from E866. The slope and absolute 
height of the curve agree with the data quite well. However, 
since here the experimental error bars are rather large, we will 
make comparisons to M spectra from other experiments in 
Secs. lIV Dl and lTVFl to test our model further. 



In this section we compare our results with the measured 
M spectrum from E866, which are basically the pj spectra 
integrated over pj. The double-differential cross section is 
given by the E866 collaboration as 



M 



da 



AM Ax i 



(132) 



Again the data are given in several bins of M and Xf and for 
every datapoint the average values (M) and (xf ) are provided. 
For the different contributions in our model we calculate the 
quantity of Eq. (1 1 32b by integrating over p 2 for every data- 
point using these averaged values: 
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(133) 



The maximal possible pj is determined by the kinematics to 



FIG. 20. M spectrum obtained from our full model. Everywhere 
D = 0.45 GeV, T = 0.2 GeV, A = 5 MeV, K q = 1 and K g = 2. The 
PDFs are the MSTW2008LO68cl set. Data are from E866 binned 
with -0.05 < x F < 0.05. Only statistical errors shown. 



P 1 + P 2 = q + X 
(Pi +P 2 ~ qf =X 2 = M 2 R 
> S + M 2 - M 2 R = 2 {P x + P 2 ) q 
= 2 V5£ 



(134) 
(135) 



2 V5 tJm 2 + (p r )Lx + 4(9*) 



(136) 



C. E772 



Experiment E772 B55I1 measured dimuon production in pd 
collisions at S ~ 1500 GeV 2 . For the calculation of the triple 
differential cross section we again use Eq. d76i > and for the 
average values of M and x F we use the center of the M and x F 
bins. Since the experiment was done on deuterium we have 
calculated pp and pn cross sections and then averaged. 
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FIG. 21. p T spectrum obtained from our full model with different 
PDF sets. Everywhere D = 0.45 GeV, T = 0.2 GeV, A = 5 MeV, 
K q = 1 and K g = 2. Data are from E772 binned with 5 GeV < M < 6 
GeV, 0. 1 < x F < 0.3. Only statistical errors are shown. 



In Figs.l2Tlandl22lwe compare the results of our full model 
with different PDF sets to triple differential data from E772 in 
different M bins. Agreement is again quite well, however, the 
shape of the spectrum seems to favor a slightly smaller value 
for D, which would enhance the spectrum near pj = 0. Never- 
theless, we have chosen D = 0.45 GeV, since this value allows 
us to describe the data from several different experiments with 
only minor deviations. 



on copper we have calculated pp and pn cross sections and 
then averaged (29 protons and 34 neutrons). 
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FIG. 23. pr spectrum obtained from our full model with different 
PDF sets. Everywhere D = 0.45 GeV, T = 0.2 GeV, A = 5 MeV, 
K q = 1 and K g = 2. Data are from E605 binned with 7 GeV < M < 8 
GeV, Xp = 0.1. Only statistical errors are shown. 



In Fig. [23] we compare the results of our full model with 
different PDF sets to triple differential data from E605. Here 
the shape of the spectrum confirms our chosen value for D and 
the overall agreement is good. 



> 

CD 

a 



LU 



0.1 



0.01 



MSTW2008lo68cl 
GJR08IO 



0.5 1 1.5 2 2.5 3 
P T [GeV] 



FIG. 22. pj spectrum obtained from our full model with different 
PDF sets. Everywhere D = 0.45 GeV, T = 0.2 GeV, A = 5 MeV, 
K q = 1 and Kg = 2. Data are from E772 binned with 7 GeV < M < 8 
GeV, 0.1 < Xf < 0.3. Only statistical errors are shown. 
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FIG. 24. M spectrum obtained from our full model. Everywhere 
D = 0.45 GeV, T = 0.2 GeV, A = 5 MeV, K q = 1 and K g = 2. Data are 
from E605 with x F = 0.125. Only statistical errors shown. 

In Fig. [24] we plot our result with different PDF sets for the 
double differential cross section together with the data from 
E605. Again agreement is quite good over the entire range of 
M. 



Experiment E605 ll56ll measured dimuon production in pCu 
collisions at S ~ 1500 GeV 2 . For the calculation of the triple 
differential cross section we again use Eq. ( T76b and for the 
average value of M we use the center of the M bin. For the 
pj spectrum E605 gives Xf =0.1. For the double differential 
cross section we use Eq. (11331 l. Since the experiment was done 



E. E288 

Experiment E288 lf57ll measured dimuon production in pA 
collisions at S ~ 750 GeV 2 . For the calculation of the triple 
differential cross section we again use Eq. ( |76*| | and for the 
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average value of M we use the center of the M bin. For the 
p T spectrum E288 gives for the rapidity y = 0.03 and we thus 
chose Xf = for our calculation. The experiment was done on 
different nuclei and only data averaged over the results from 
these nuclei have been presented. Therefore, we have calcu- 
lated pp cross sections only. 



our results for the double differential cross section with the 
data from E439. For both PDF sets the absolute height and 
slope agrees well with the data. 
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G. E537 (Antiprotons) 



Experiment E537 ll58ll measured dimuon production in pW 
collisions at S ~ 236 GeV 2 in an invariant mass range of 
4 < M < 9 GeV. The obtained cross sections are double dif- 
ferential in two of the observables M, x F and p\. To calculate 
the cross sections differential in p 2 T with our model we use 



FIG. 25. pr spectrum obtained from our full model with different 
PDF sets. Everywhere D = 0.45 GeV, T = 0.2 GeV, A = 5 MeV, 
K q = 1 and K g = 2. Data are from E288 binned with 7 GeV < M < 8 
GeV, y = 0.03, we have chosen x F = in our calculation. Only 
statistical errors are shown. 



da 



dxfdpj 



I 

Jm 2 



d(T 



bin dM 2 dx F dpl 



-dM l 



da 



dM 2 dx F dp\ 



((M i ),(x F ),(PT)) ■ (140) 



In Fig. [25] we compare the results of our full model with 
different PDF sets to the triple differential data from E288. 
Again the agreement with the data is good and confirms our 
choice of parameters. 



F. E439 



The sum runs over several mass bins, which we choose as 
M = 4 ... 5, 5 ... 6, 6 ... 7, 7 ... 8, 8 ... 9 GeV and in each bin 
we take the central value for {Mi). Since the experiment was 
done on tungsten we have calculated pp and pn cross sections 
and then averaged (74 protons and 1 10 neutrons). 



> 

CD 

jo 



"D 

"6 

"D 



0.1 



0.01 



0.001 



MSTW2008lo68cl 
GJR08IO 



4.5 



5.5 6 6.5 
M [GeV] 



7 7.5 



FIG. 26. M spectrum obtained from our full model with different 
PDF sets. Everywhere D = 0.45 GeV, T = 0.2 GeV, A = 5 MeV, 
K q = 1 and Kg = 2. Data are from E439 with x' F = 0.1. Only 
statistical errors are shown. 



The details of experiment E439 and how we calculate the 
cross section are given in Sec. HIE 21 In Fig. [26] we compare 



We compare the calculated pr spectrum with the data in 
Fig. [27] Our full model is on the lower side of the error bars 
of the data. However, one should note that the experimental 
error bars are rather large and thus the possibility to confirm 
or rule out our model is limited. 

To calculate the M spectra we use 



dcr 



dMdx, 



-f 

Jo 

-f 

Jo 



dp 2 T 



da 



dp z T 2M 



dMdx F dpj 
da 



dM 2 dx F dp 2 



(M,x F ) , (141) 



with (pr)max gi ven m Eq- d 1 39b . We compare our calculated 
M spectrum with the data in Fig. [28] Agreement is better than 
for the pj spectrum, however, the experimental error bars are 
again large compared to proton data. 
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FIG. 27. p T spectrum obtained from our full model. Everywhere 
D = 0.45 GeV, T = 0.2 GeV, A = 5 MeV, K q = 1 and k s = 2. The 
PDFs are the MSTW2008LO68cl set. Data are from E537 with 4 
GeV < M < 9 GeV, < x F < 0.1. We have chosen x F = 0.05 in our 
calculation. Only statistical errors are shown. 



use a modified version of Eq. ( l76l >: 



2 VS£ 



dcr 



2 



tt(S - M 2 ) dx^dp 

2 VS£ 
^ ;r(S - M 2 ) 

2yfSE 
~ tt(S - <M> 2 ) 



dcr 



bin dM 2 dx' F dpj 



dM z 



AM 



do- 



dM 2 d^dp: 



■{{M),x' F ,p T ) , (142) 



where 



£= ^{M) 2 +p 2 T +{x' F ) 2 {(q'^) 



(143) 
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FIG. 28. M spectrum obtained from our full model. Everywhere 
D = 0.45 GeV, T = 0.2 GeV, A = 5 MeV, K q = 1 and Kg = 2. 
The PDFs are the MSTW2008LO68cl set. Data are from E537 with 
< Xp < 0.1. We have chosen x F = 0.05 in our calculation. Only 
statistical errors are shown. 



and AM 2 = M 2 m „ - M^ in with M max (M^) the upper (lower) 
limit of the bin. For the average value of M we use the center 
of the M bin and we choose everywhere x' F = 0. In Figs. l29ll30l 
we show our predictions for different values of F in different 
M bins. Note that while at E866 energies we could not dis- 
criminate between different T over a wide range (cf. Fig.fT2l). 
the results here become more sensitive to this parameter. 
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FIG. 29. pr spectrum obtained from our full model for different 
values of T. Everywhere D = 0.45 GeV, A = 5 MeV, K q = 1 and 



K g = 2. The PDFs are the MSTW2008LO68cl set. 
GeV < M < 2.5 GeV. 
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H. Prediction for PANDA 



Based on the parameters which we have fixed on the avail- 
able data above, we here present our predictions for DY pair 
production at S =30 GeV 2 in pp collisions, where, for exam- 
ple, PANDA |0] will measure. 

For the calculation of the triple differential cross section we 



In Fig. [5T| we compare results for different PDFs and the 
uncertainty induced by the choice of different PDFs is com- 
parable to the uncertainties we found for high energies (e.g. at 
E772). To ensure that at such low hadronic energies the ficti- 
tious gluon mass is still small enough, we study our model at 
different values of A, see Fig. [32] The results coincide, indicat- 
ing that our standard choice of A = 5 MeV is still applicable 
at these energies. 
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FIG. 30. spectrum obtained from our full model for different 
values of T. Everywhere D = 0.45 GeV, A = 5 MeV, K q = 1 and 
K g = 2. The PDFs are the MSTW2008LO68cl set. 4 = and 2.5 
GeV < M < 3.5 GeV. 



at the low hadronic energies of the PANDA kinematics, we 
again vary one of the two parameters and keep the other one 
fixed and show our results in Figs. [33] and [34] The results 
for different K q deviate by about 15%, which is comparable 
to the deviation at E866 energies. However, the results are 
practically insensitive to variations in K g . 
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FIG. 31. pr spectrum obtained from our full model for different PDF 
sets. Everywhere D = 0.45 GeV, T = 0.2 GeV, A = 5 MeV, K q = 1, 
K . = 2. x' = and 1.5 GeV < M < 2.5 GeV. 
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FIG. 32. p T spectrum obtained from our full model for different 
values of A. Everywhere D = 0.45 GeV, T = 0.2 GeV, K q = 1 and 
K g = 2. The PDFs are the MSTW2008LO68cl set. 4 = and 1.5 
GeV < M < 2.5 GeV. Note that both curves are basically on top of 
each other. 

To check the dependence of our results on the choice of the 
subtraction parameters K q and K g (see Sec. IIIIEI for details) 



FIG. 33. pt spectrum obtained from our full model with differ- 
ent values of the subtraction parameter K q . Everywhere D = 0.45 



GeV, T = 0.2 GeV, A = 5 MeV and k s 



2. The PDFs are the 



MSTW2008LO68cl set. x' = and 1.5 GeV < M < 2.5 GeV. 
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FIG. 34. p T spectrum obtained from our full model with differ- 
ent values of the subtraction parameter K g . Everywhere D = 0.45 



GeV, T = 0.2 GeV, A ■■ 
MSTW2008LO68cl set. 



5 MeV and k„ 



1. The PDFs are the 



and 1.5 GeV < M < 2.5 GeV. Note 



that the curves are practically on top of each other. 

Finally in Fig. [35] we compare our predictions with 
a PYTHIA calculation (PYTHIA version 6.225, CTEQ5L 
PDFs), each for different values of the average initial kj. As 
explained above, PYTHIA calculations for E866 conditions 
seem to prefer a somewhat smaller width for the initial kj dis- 
tribution, (k 2 T ) = (0.8 GeV) 2 instead of (0.9 GeV) 2 . Since var- 
ious calculations ([59], [60]) hint to some monotonic depen- 
dence of the initial kj with the underlying VS~ , we would ex- 
pect that at PANDA energies, a somewhat smaller (kj) should 
be used. Therefore, in Fig. [35] we also show calculations 
with (k 2 T ) = (0.6 GeV) 2 for PYTHIA and (k 2 T ) = (0.7 GeV) 2 
(D = 0.35 GeV) for our model. For PYTHIA, already with 
(kj) = (0.6 GeV) 2 , the difference in the functional behaviour 
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is rather large compared to the PYTHIA calculations with the 
higher parameter values. This may be taken as some hint for 
the theoretical uncertainties. On the other hand, the intrinsic 
kj in PYTHIA is some effective parameter. It is not clear, 
whether this parameter should follow the same energy depen- 
dence in pp and in pp collisions, since multiple effects are en- 
coded. Note, that the PYTHIA results shown in Fig. [35] were 
not multiplied by a K factor, i.e. K = 1. 
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FIG. 35. pr spectrum obtained from our full model and from 
PYTHIA (see main text for details) for different values of the av- 
erage initial k T . Our predictions were calculated with T = 0.2 GeV, 
A = 5 MeV, K q = 1, Kg = 2 and the MSTW2008LO68cl PDF set. 
x' = and 1.5 GeV < M < 2.5 GeV. 



V. CONCLUSIONS 

In this paper we have extended a phenomenological model 
of DY pair production We have included all relevant 

processes up to NLO in the strong coupling a s to account for 
the missing strength in the LO calculation (K factor). To de- 
scribe DY pair pr spectra we introduced initial parton trans- 
verse momentum distributions and to regularize the otherwise 
divergent pj spectra at NLO we distributed the masses of the 
quarks with spectral functions. To avoid double-counting we 
introduced a subtraction scheme which removes those 0(a s ) 
contributions which are usually absorbed into the renormal- 



ized quark PDFs. 

The results show that with our choice for the width of the 
initial parton transverse momentum distribution, D as 0.45 
GeV ((kj) as (0.9 GeV) 2 ), the shape of the pj spectra is repro- 
duced very well. We note, however, that this width might be 
S dependent, which introduces additional uncertainties. The 
spectral functions (quark mass distributions) serve as a regu- 
lator for the NLO order processes, which for massless quarks 
are divergent for pj — > 0. However, at high energies (for 
example at E866) over a wide range the results depend only 
weakly on the width F of the mass distribution. In addition we 
could fix the subtraction parameters K q and K g at E866 energies 
and found that they are of natural magnitude (0(1)). 

In summary we found that in our phenomenological model 
we can reproduce measured pj and M spectra for DY pair 
production in pp and pA reactions at different hadronic en- 
ergies without the need for a K factor. The comparison to 
pA data from E537 is inconclusive, however, our full model 
is still within the large (compared to the pp and pA data) er- 
ror bars. The general agreement of our results with the data 
from different experiments indicates that we effectively have 
parametrized the soft initial state interactions in the nucleon 
by fixing our parameters on the available data. Using this 
framework we have obtained predictions for DY pair produc- 
tion at the low hadronic energy regime, where future experi- 
ments, for example PANDA, are aiming at. We have found 
that our predictions become more sensitive to the mass distri- 
bution width F, which we could not reliably fix at higher ener- 
gies. In addition we found some sensitivity on the subtraction 
parameter K q , which is comparable to the finding at high ener- 
gies (E866). Nevertheless, our model provides a narrow band 
of estimates for the fully differential DY pair production cross 
section at low energies. 
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Appendix A: Gauge invariance 
1. Electromagnetic sector 

By assigning to the annihilating quark and antiquark different masses m\ and m 2 , see Fig. [3 gauge invariance is broken at the 
quark-photon vertex. One can easily see this by contracting the quark current with the photon momentum q^. 

v(p 2 ,m 2 )y u(p u mi) ■ 
= v(p2,m 2 )qu(pi,m]) 
= v(p2,m 2 )(f l + f 2 )u(p u mi) 

= v(p 2 ,m 2 )(m l - m 2 )u(pi,m l ) + . (Al) 
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However, in the full amplitudes for DY pair production of Secs.[II]and[lII]this is not an issue since gauge invariance is preserved 
at the lepton-photon vertex. To realize this one has to look at the gauge dependent part of the photon propagator. The propagator 
has the following Lorentz structure: 

G(?)~(^ y -^fj, (A2) 

with a gauge parameter ^. Now we insert the gauge dependent part of Eq. ( 1A2I ) between the quark and lepton currents and exploit 
the Dirac equation: 

v(p2,m2)y"u(pumi) ■ (q^q v ) ■ u{k\,m)y v v{ki,m) 

- v(p2, mi) f u(pi, mi) ■ u(k\,m) $v(k.2,m) 

= v(p2, mi) (p x + f 2 ) u(pumi) ■ u(ki,m) (ft + fc) v(fe, m) 

— v(f2, mi) (mi — mi) u{p\,m\) ■ u(k\,m) (m — m) v(fe, tri) 

= , (A3) 
and so the amplitude is invariant under gauge transformations of the electromagnetic field. 



2. Strong sector 



As shown in Figs. [6] [7] [8] and [9] we always keep the quark masses fixed at any quark-gluon vertex. This ensures that all 
our calculations are also gauge invariant in the strong sector. This is trivial for the gluon bremsstrahlung and gluon Compton 
scattering processes. For the vertex correction (VC) and self energy diagrams the proof proceeds as follows: We denote the LO 
(a° s ) DY amplitude with Mlo and the vertex correction amplitude with Myc- Then 

M h0 ~ v(p2,mi)f u(p\,m\) ■ L M (q 2 ) , 



M< 



v(p2,m 2 )r f '(q )u(p u m l ) ■ L^q ) 



(A4) 
(A5) 

where is the leptonic part and V is of order a s . Also of order a s are the interferences of the LO process with the self energy 
diagrams in Fig. |7]and we have to include them by dressing the external quark legs with field strength renormalization factors of 
VZ2(m,) with Z2 = 1 + 6Z2 + 0{a 2 s ). Then to order a s the amplitude for the process qq — > Z + T can be written as 

M = yjZ2(m x )ylZ2{m2)-v(p2,m2) + r"(q 2 )) u{p u m x ) ■ L^q 2 ) 

\ + l -5Z 2 (m x ) + 0(a£)J|l + l -5Z 2 (m2) + 0(a*)J • v(p 2 ,m 2 ) (xe q y" + T"(q 2 )) u{p u m{) ■ L^q 2 ) 

ie q y |l + l -6Z2{ mi ) + i« 2 (m 2 )j + rV) 



= v(p 2 , m 2 ) 



u(p\,m{) ■ Lfiiq ) + 0(a s ) 



(A6) 



To prove gauge invariance we now have to calculate the gluon gauge dependent parts of and 6Z2, which we denote by the 
index g. We begin with Pi. We denote the gluon momentum with k, insert only the gauge dependent part of the gluon propagator 
(k°kP) and again exploit the Dirac equation: 

-14 



i(-p> 2 -f + m 2 ), i 



v(p 2 ,mi)Ptu(pum{)= \ d ei kv(p 2 ,m 2 )(igy a f) 

J (P2 + k) 2 - m\ 

- #, - ft + m2 
G> 2 ,m 2 )#— *' 



)u(pi,m 1 ) ■ —K2~T2 



eAna s (t 



"t") J d A kv{ 



(pi - k) 2 - m\ 

j - ft + m\ 



2 T ~ ' — 7 2 fu(p umi )- — 

(p 2 + k) 2 -m\ (pi- k) 2 - m\ 



1 

F 



= - e g 4na s (ff) \ d 4 k v(p 2 , mi) ii-f 2 - mi) - {-p 2 - ft - mi)) 

J (P2 + 



f 2 -f + m 2 



— t + mi 



■f- 

(P2 + k) 2 ~ ml {pi - k) 2 - m 2 



x u(p u mi) ■ 

/§, — f + m\ , > 1 

d 4 kv(p2,m 2 )(-l)y j r +f + m l )-(-p 1 +OTi)j u(pi,mi)--r 

(pi - k) L - OTj ' ' k? 



e a 4na s (ft 



e q 4na 



a ) J d 4 kv> 
Jiff) J d 4 kv> 



(p2,m2)(-l)y(-l)w(pi,mi) ■ 

y 

(p2,m 2 ) -jr u (Pu m i) ■ 



1 



(A7) 
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The renormalization factor 5Z\(m) is connected to the QCD quark selfenergy via [43 



6ZUm) 



dW(p) 



The gauge dependent part of the quark selfenergy is given by 

-tV(j>) = J< : 



(pi - k) 2 - m 2 



ik a kP 1 

k 2 ¥ 



/. p,-t + mi 1 
d 4 kf-/ ' ■ 



"(pi-£) 2 -m, k4 



Now we find for the renormalization factor: 
6Z%(m)=i4na s (ft a ) j d 4 kf- 



■ iAna s (t 



d 4 kt- 



= iAnaJt 



— iAna 



Thus 



: iAna s (t 



v(pi,m 2 ) 



a f) J 

a f) J d 4 

s {ff) j d 4 
a f) J d 4 k-^- 



(p — k) 2 — m 2 
1 



f=m 



k 4 



-t + m){-2f + It) 



- kip + It 2 - m 2 



'> 2 _ m 2 



k 4 



1 (-£ + 2m)(-2m + If) 

+ 



■2mt + t (_ 2m f + f} 

' -2mt + k 2 + (-* + 2m)(-2m + 2f) 
k 2 (Am 2 - Ami + k 2 ) 2 



1 

' k 2 
1 

k 2 



ie q y -5Z\{m x ) + -6Zl(m 2 ) + F 



ie„y iAna s 



u(pi,mi) 



d 4 k 



-1 



= v(p 2 , m 2 ) 
= 

and so the amplitude in Eq. ( IA6b does not depend on the gluon gauge 



k 4 ' q 



siff) j 



4,7" 



d«k 



k 4 



u(p u mi) 



(A8) 



(A9) 



(A10) 



(All) 



Appendix B: F t 

As already mentioned in Sec. IIII A 1 1 we unintentionally renormalized the charge at the quark-photon vertex by assigning 
different masses m\ and m 2 to the annihilating quark and antiquark, which breaks gauge invariance at the quark-photon vertex, 
see Appendix IA II Thus 

\m\F x (q 2 ,m\,m\) + 1 , (Bl) 

which we illustrate in Fig. [36] There we plot the real part of the correction 6F\ to F\ to order a s for small -sjq 2 . The correction 
is defined by 

Fi = 1 + ^ SFt . (B2) 
An 

As one can see, 6F\ approaches zero for the case of equal quark masses nt\ — m 2 — 0.1 GeV, as it should. However, for different 
quark masses (in our example plot: m\ =0.1 GeV, m 2 = 0.5 GeV) this is clearly not the case. 

This behavior could potentially spoil our calculations. One should note, however, that we always stay far away from q 2 = 0, 
since q 2 = M 2 sets the hard scale for our calculations. Therefore, reasonable physical values of q 2 are larger than 1 GeV. To 
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FIG. 36. Correction to F\ at order a s for equal quark masses (solid) and different quark masses (dashed). See main text for details. 



study the influence of different quark masses in the physically interesting range of q 2 we devise the following scheme: 
to calculate the hadronic cross section we weight the partonic subprocess cross sections by quark mass distributions (spectral 
functions), see Eqs. (l71H73t . Thus also the form factor F\(q 2 ,m 2 ,m^) is weighted by these distributions in our calculation. 
Therefore, it is worthwhile to compare the weighted form factors for different masses, F\{q 2 , m 2 , m 2 ), and for equal masses, 
F\ (q 2 , m 2 , m 2 ), for physically interesting q 2 . Because of Eq. (IB 2b it suffices to compare only the corrections 5F\ . We define 

r m l 

6Fi(q 2 )= dm 2 A(p)6Fi(q 2 ,m 2 ,m 2 ) , (B3) 
Jo 

with the spectral function A(p) defined in Eq. ([72}. Now we know that 8F\(q 2 , m 2 , m 2 ) shows the correct low q 2 behavior, 

lim 6Fi (q 2 , m 2 , m 2 ) = , (B4) 

and we know that the spectral function is normalized to 1, see Eq. d74l >. Therefore, also the weighted correction shows the right 
behavior for q 2 — > 0: 

r m N 

lim 6F x (q 2 ) = dm 2 A(p) ■ lim 6Fi(q 2 , m 2 , m 2 ) = 1 • = . (B5) 

f^o Jo q^0 

Next we define the weighted correction for different masses: 

r m l „ r m l „ , . „ 

5Fi(q 2 )= dm\ dm 2 1 A(p 1 )A(p 2 )6F l (q 2 ,m 2 l ,m 2 2 ) . (B6) 
Jo Jo 

In Fig. [37] we compare the real parts of 5F\{q 2 ) and 6F\{q 2 ) for yjq* =1 ... 20 GeV. As one can see they agree very well 
over the entire range. Thus we conclude that the wrong behavior of SF\(q 2 , m 2 , m^) near q 2 = ultimately does not affect our 
calculations. 



Appendix C: Kinematics 

The calculation of the hadronic cross sections for the LO process, see Sec.[IIJ gluon bremsstrahlung, see Sec. lIIIBl and gluon 
Compton scattering, see Se c. IIH CI requires to remove unphysical solutions for the momentum fractions x,. We closely follow 
the arguments presented in 027ll and refer to this publication for more details. 



m 1 =m 2 =0.1 GeV 
m 1= 0.1 GeV, m ? =0.5 GeV 
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FIG. 37. Comparison of the real parts of the weighted corrections to F\ at order a s for equal quark masses (solid) and different quark masses 
(dashed). For the spectral functions a width T = 0.2 GeV and a large quark momentum component p + = | was chosen, cf. Eqs. d7 1 173t . The 
gluon mass was set to A = 0.05 GeV, cf. Eq. J95b . Both curves agree well over the entire range. 



1. Bremsstrahlung 

We begin with the symmetric case of gluon bremsstrahlung. The hadronic cross section reads, cf. Eqs. ( 12315 It , 
dcr B 



dM 2 dp 2 dx F 



■ -j~ dxi -j~ dx2 j dpi x J" dp2 ± J~ dm 2 J" dm 



^ qffi{xi,Pi ± ,m\, q 2 )fi(x 2 , p2 ± ,m\, q 2 ) 



2 Vsp cm (q z ) mdx do- B 



dM 2 df 



(CI) 



■6(( Pl +p 2 -q) 2 -A 2 ) . 

Again, the f { are our unintegrated parton distributions, see Sec. lIID2l the partonic cross section d ^" 2 B d) is given in Eq. ( 110U and 
A is the fictitious gluon mass, introduced to regulate the soft divergence. Now we collect everything except 5- and ©-functions 
in F and rewrite the cross section as 



dcr 



dM 2 dx 



*—— = \ dx\ \ dx 2 I dp l± dp 2± dm\ \ dm 2 2 F(x u p\^m 2 ,x 2 ,p2 ± ,m 2 2 ,M i 
Fdpj Jo Jo J J J J 

x 6 (( Pl + P2 -q) 2 -A 2 ) &(E g ). 

I 



(C2) 



The 5-functions in Eq. (IC2t must be worked out in a way that Adding and subtracting Eqs. (IC5b and (IC6b yields 
allows to discern physical and unphysical solutions for the 

momentum fractions x, in order to perform the ^-integrations. ^ 
For this aim it is useful to rewrite the parton momenta in terms 
of different variables: 



1 m\+m\ 

-—q H 

A H 2 



q = pi + p2 , 



(C3) 
(C4) 



k ■ q = 

H 2 

Solving Eq. dC7| > for Ic + yields 



Inverting the last two equations, we can use the on-shell con- 
ditions for the partons to obtain 



P J 2 m\+m\ 



(C7) 
(C8) 

(C9) 



1 



1 



m 1 =p 1 =\—q — k\ — — q —k-q + k 



and 



Inserting this result into Eq. dC8b leads to an equation 
'■ 2 (C5) quadratic in k~ : 

m 2 — m 2 — k + q~ + kTq + — 2k ± ■ (q ± ) 



1 



1 



ki-\q 2 + '± 



m 2=P2 = \^ + k \ = 7$ +k-q + k 



(C6) 



-q- + k-q + - 2k ± ■ (q ± ) (C10) 
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= (k~) q + + k~(-2k ± ■ (q ± ) -m\+ m\) 



+ \kl-79 + ~ 



(Cll) 



and 



Pi 



1 



-A 



(C15) 



The solutions are 



2§+ 



q + 2q + 



C - 1 2 



(C12) 



Since there are two solutions for k~ and k + , respectively, we 
also get two solutions for x\, x-i- To determine which set of 
x\,x<i and thus k + ,k~ has to be chosen we take the limit of 
zero parton transverse momentum and vanishing masses, 



(C16) 



(C17) 



Inserting dC12| > into dClOl l gives the solutions for k + : 



Inserting expressions ( IC16b and ( lC17b into ( IC141 > and ( IC15b 
yields two solutions for the momentum fractions, 



1 



k± ■ (q±) , m 2 



r 



2q+ 



+ q - 



4 q kl 2 



(C13) 



Rewriting Eqs. ( f27l > and (f28j) in terms of q and k we obtain the 
solutions for the parton momentum fractions (P 2 = P\)'- 



(C14) 



and 



if 



(C18) 



(C19) 



The upper solutions correspond to the unphysical case x\ = 
X2 = 0. Thus we only keep the lower solutions when evalu- 
ating the phase space integrals. This requires the integrals in 
Eq. ( lC2b to be evaluated in the correct order, otherwise one 
cannot disentangle the two different solutions for x\ and X2- 

We begin by introducing several integrals over ^-functions 
in Eq. flC2b . In this way we will transform the integration 
variables to the above chosen q and k±: 



dM^dxdp 2 = jo dX ' f I I ^ 2± I d ^ ± ^J" dm ' f dm l J d <7 + J d <7 F(x u pi ± ,m 2 ,x 2 ,p2 ± ,m 2 ,,M 2 ) 

XS(q + - (p\ +P + 2 ))6 (§~ - (Pi + P 2 )) ^ " (Pi x + P2 ± )) 5 (2) - \ (p l± - frjj 

x 5 ((pi + p 2 - qf - A 2 ) (E g ) . (C20) 



First we perform 



j dp ix j Ap 2 J {1) ((§7) - (p l± + P2j) 5 (2) (4 - \ (P2 ± - Pi j) = 1 ■ (C21) 



Now we calculate the integral 

T d^i t dx 2 5 (q + - {p\ + p 2 )) 5 (q~ - {p\ + p 2 )) . 
Jo Jo 

(C22) 

According to Eqs. ( IC 1 2b -( IC 15b the ^-functions in this expres- 
sion have two possible solutions for each p\ and p 2 . How- 
ever, as explained above, we now have to explicitly remove 
the unphysical solutions (xi)+ and (xz)- , which are the ones 
corresponding to the upper sign in Eqs. dC 12b and dCl 3b : 
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■f dx, -f dx 2 S(q + -{p\ + p 2 ))6(q -(p l + p 2 )) 
Jo Jo 



Jo Jo 



= r in r 

Jo Jo 



dx 2 5 



dx 2 <5 



(i(§J-4) 2 



+ m, 



+ mf 



(n) 2 



+ to; 



x 2 P- 



q -X\P 1 



2 2 ^ 



q -(xi)_Pj 



2 \ 
+ m 2 



(^2) +J P7 



(^i)'(^(pr) 



■ 0(1 - (*0_) 0((*i)_) 0(1 - (x 2 ) + ) 0((x 2 ) + ) . (C23) 



Using dq + dq = 2dqodq z we can evaluate one of the remaining integrals of Eq. ( lC20b with the help of the <S-function: 



2 J dq Q S(( Pl +p 2 -q) 2 -A 2 ) = 2 J dq 6((q-q) 2 -A 2 ) = -^ 



(C24) 



with Eg = VO 3 ) 2 + A 2 = *\j((q) -4} Collecting the pieces, Eq. ( IC20I ) simplifies 



to 



do- 



dM 2 dx 



dk ± I dm 2 



= d M m,x d ^) 

fdpj, J(q^ n JO JO 



dm 2 F((JC!)_, to 2 , (x 2 )+, /? 2± , to 2 , M 2 ) 





{P- X f- 



(x0 2 (x 2 )i(p r ) 



x 0(1 - (xi)_) 0((xi)_) 0(1 - (* 2 )+) ® ((**)+) . 



(C25) 



Now (xi)_,p I± , (x 2 )+ and /? 2± are fixed: 



fe)+ 



2 ™2 



r fe± ■ (gx) ot 2 - 

2 2§+ + \ 



k ± -(q ± ) i m 2 - 



2 m 2\ 



r 



-q 2 -ki 
\ H ± 2 



1 



q + k ± ■ (q ± ) m 2 - 



2 ™2 



+ 

2 § 



24- \ 



£l-(£J W 2~ W 1 ^ 

§ 2§" 



fi 2 g2 ^ + ^ ) 



P2, = ^(fj + £l , 

k± ■ (q±) = k ± q ± cos <f>k ± 



(C26) 

(C27) 

(C28) 
(C29) 
(C30) 
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with 



q + — Eg + E g + q z , 
q~ =E q +E g - q z , 



E q = tJm 2 +p 2 T + q\ , 

Eg = a/(§.l) 2 + q z - 2(?±) • q± - 2q~ z ■ qz + p 2 T + ql + X 2 . 



(<7.l) • q± = qxpr cos <f>$ x , 



(C31) 
(C32) 

(C33) 

(C34) 

(C35) 
(C36) 



The integration limits can now be found from general considerations: (m^max is fixed by the condition that and (^2)+ must 
be real numbers, 



(WjW = -2&l • (q±) + m\ + q + q - \\4q + q m\ + 4q + q | ^(q±) - k ± 



From (m^max > we find 



OhW = 2k ± ■ (q ± ) + q q - 2 \\q + q I — + &o 



and from (m 2 ) ma x > follows 



1*0. 



2 

max 



q + q (q + q -q 2 ±) 



4 q\ cos 2 (<f> k j) 

The energy of the incoming partons must be less than the energy of the hadronic system, qo < VS~, and so 

I (lo3 La = Pt cos 4>q + J p\ (cos 2 (pq - l) - (q z - q z ) 2 + ( Vs - E q f - A 2 . 
Finally, q ± is a real number and thus 

tizG = Qz ± ^P 2 T {cos 2 ^-\) + {<S-E c ) 2 -A 2 . 



(C37) 



(C38) 



(C39) 



(C40) 



(C41) 



2. Leading order process 

The kinematics of the LO cross section is a special case of the bremsstrahlung kinematics of Appendix IC 1 1 Namely at LO 
the four-momentum of the incoming partons is equal to the four-momentum of the virtual photon: q-q. From Eq. ( f52l i we note 
for the partonic cross section 



d<r 



LO 



5 (M 2 - (p, + P2 f) 6 [p\ - {p l± + p 2 J 2 ) 6 ix F 



(Pl)z + (P2% 



(C42) 



dM 2 dx F dp 2 T ' v ~ ~ ' r "' 1 ' v 1 ~ ' x "' x ' ' ' \" (q z ) n 
Now employing Eqs. JC4HC19b and Eqs. dC2UlC23l and everywhere replacing q by q we find for the LO hadronic cross section 
dcr LO 



dM 2 dx F dpl 



J 2dq {) J dq z J dq ± J dlc ± J dm 2 j dmlF LO ((xi)-,pi ± ,m 2 ,(x2)+,P2 ± ,ml,M 2 ) 



{\q±-k±f +m\ (±<? x +*Yf + 



X (1 - ((*)_) 0(1- (x 2 ) + ) ((x 2 )+) 



:6(M 2 -q 2 )6(p 2 T -(4 ± ) 2 )6\x F 



q z 

(<7z)ma; 



(C43) 
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With help of the three remaining 5-functions we can now easily perform the four integrations over the components of q: 



Act 



r-' T r$l'%mx 1 {W\Kmi mfiu w(n\ 

= f d 0± f IdO?J» f dm? f dm* 2^S 
fdpt. Jo Jo 2 Jo Jo Eg 



AM 2 Ax 



( p i) 



[\q±-k^f +m\ [\q,+k^f + 



(Xl)l(x 2 f + (Pi) 

x0d-(n)-) ®((*i)-) ®d -(x 2 ) + ) ®((x 2 ) + ) . 
The integration limits are now recovered from Eqs. ( IC37HC39b , again by replacing q with q. 



^lo(Ui)-, p\ x ,m\, (x 2 )+, P2 ± ,mj, M 2 ) 



(C44) 



3. Compton scattering 



At this point we like to stress the kinematical differences between bremsstrahlung and Compton scattering. For bremsstrahlung 
we have a quark from nucleon 1 annihilating with an antiquark from nucleon 2 or vice versa. However, we treat quarks and 
antiquarks on equal footing and distribute their masses with the same spectral function, cf. Sec. Ill D 21 Thus we can easily take 
care of both cases by simply summing over all quark- and antiquark-fiavors in Eq. ( ICU . Gluon Compton scattering is different 
since we keep the gluons massless and the simplification from above does not apply anymore. However, we can calculate one of 
the two cases, for example quark/antiquark from nucleon 1 annihilates with gluon from nucleon 2, and then simply find the other 
case by symmetry considerations: nucleon 1 and 2 are defined by their direction of motion along the z-axis. Thus by changing 
z to -z and so xp to —Xp we find that the second case corresponds to the first case with xp — > — Xp. The hadronic cross section 
therefore reads, compare with Eq. ( IClt , 



AM^Tp = { ±Xl f ^ 1^ f** 



qtifdxixupi^m^q )g 2 (x 2 ,P2 ± ,q ) - AM 2 At ^ Pl 



qf) 



J dxi j- Ax 2 J Api^ J Ap 2x J" Am 



V 2/$\ I -> 2 2\~ , -* 2\ 2 V*Pcm(?z)max d(Tc 



x 2_i<li{jd2{xupi x ,m v q )gi(x 2 ,p2 ± ,q ) 



dM 2 dr 



5 ((pi +P2- qf) 



(dcr c )i2 (dcr c )21 



AM 2 Ap 2 Ax F AM 2 Ap 2 Ax F 



(C45) 



The indices 1 and 2 for the parton distributions denote the parent nucleons (p,n,p). g is the transverse momentum dependent 
gluon distribution function and we choose it in analogy with the transverse momentum dependent quark distribution function of 
Eq. d67]), 



g(x,p ± ,q 2 ) = g(xi,q 2 ) ■ f±(p i± ) 



(C46) 



with f ± defined in Eq. d68l > and with the usual gluon PDF g. Now can we proceed similarly to Sec. IC II 



(do-c)i 



i}M\\ k Ax~ ~ f dXl jo ^ f dPl± f dP2± f f d ^ f dm i f d <l + f d § F(x u p\ ± ,m 2 ,x 2 ,p 2± ,M 2 
X6(q + - (p\ + p + 2 )) 6(q- - (p\ + p 2 )) 6 (2) ((§3 " (Pl, + Plj) 5 {2) [t ± - i (p l± - pA 
Xd^pi + p 2 -qf-m 2 ) . 
We use Eqs. dC2Tb and (IC231 and 

J Am] 6 ((pi + p 2 - qf - m\) = 1 



(C47) 



(C48) 
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to find 



(dcrc) 



dM 2 dp 2 



2~ — = I d^ I ""d^jj | 2d<? | d£ L F(xi,^i_ L) mJ ) *2,$2 J .>Af 2 
' r dx F J( 9 ,) min Jo J(?o)mi„ Jo 



-\2 



(PT) 



(xOife) 2 ^) 

x 0(1 - (*)_) 0((xO_) 0(1 - (x 2 ) + ) 0((x 2 ) + ) 
Now (xi)-,p\ ± , (x 2 )+ and /? 2± are fixed by 



(*0- = 



q k± ■ (gj 

2 r 



(X2)+ 



1 

P7 



q + k ± ■ (q ± ) m, 

1 h -\ 

2 q~ 2q~ \ 




q- 2q 



+ 



( Iq 2 -P-^ 



Pi ± = -jiq±) ~ k± , 

P2 ± = ^(17) + k± , 
k± ■ (q±) = k ± q ± cos <f> k± , 



with 



-q±j - 


(qz - qzf 


q + 


= qo+qz, 


q 


= qo-qz, 




= yjM 2 +p 2 T +q 


(§7) • q± 


= q±PT COS (f> q± , 


qz 


= XF(q z ) mdx . 



(C49) 



(C50) 

(C51) 

(C52) 

(C53) 
(C54) 
(C55) 



(C56) 
(C57) 

(C58) 

(C59) 
(C60) 

The integration limits can now be found from general considerations. |fej.| ma x is fixed by the condition that (xi)_ and (x 2 ) + must 
be real numbers: 



\ki 



(q ± ) cos 4> k± m 



2(q + q - §2 C os 2 <f> k J ^ 



(q ± ) cos 4> k± m 



2(q + q -q 2 ± cos 2 <p k± 



7 ~ g (g g - gi) + g g 2. 
q\ cos 2 ti - q + q- 



(C61) 



From < m 2 < mjj one finds 



(<?o)min = E q + a/(^I) 2 + q\ ' ~ 2(H) • q± - 2q z ■q z +p 2 T + q 2 , 

(<?o)max 



E q + V(g^) 2 + q\ ' ~ 2(^7) • q± - 2q~ z ■ qz + p 2 T + q\ + m 2 N . 

Since the energy of the incoming partons cannot be larger than the hadronic energy, we have q < Vs and thus 
|(§j| max = Pt cos (pg ± + ^pl (cos 2 (pg ± - l) - (q z - q z f - m\ + ( V5 - Egf . 



(C62) 
(C63) 

(C64) 
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Finally, q ± is a real number and thus 



= 4z± ^P 2 T (cos 2 ^ - l) + ( Vs - E q ) 2 - m 2 N . (C65) 
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